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SUMMARY 


A  computer  program  has  been  developed  to  compute  the  inviscid  supersonic 
flow  over  delta  wings  and  smooth  wing-body  configurations.  A  second-order 
accurate  predictor -corrector  finite-difference  scheme  is  used  to  integrate  the 
three  dimensional  Euler  equations  in  regions  of  continuous  flow.  Bow  shock 
and  crossflow-induced  embedded  shocks  are  explicitly  computed  as  discontinui¬ 
ties  which  simultaneously  satisfy  the  characteristic  and  the  Rankine-Hugoniot 
conditions.  In  computing  the  flow  about  complex  wing  cross  sections,  the  use 
of  appropriate  conformal  mappings  were  an  important  factor  in  developing  a 
computational  mesh  capable  of  resolving  the  large  flow  gradients  that  are  inevi¬ 
table  in  the  vicinity  of  wing  leading  edges. 

Geometry  programs  were  developed  to  supply  the  appropriate  geometric 
boundary  conditions  necessary  to  compute  complex  wing  cross  sections.  New 
starting  solutions  were  developed  that  were  more  appropriate  to  thin-wing  cross 
sections. 

In  the  initial  stage  of  development,  the  flow  about  conical  wings  was  com¬ 
puted  for  a  two-fold  purpose.  Conical  Euler  solutions  would  be  used  as  start¬ 
ing  conditions  for  conical  wing- bodies  and  three  dimensional  wings  and  the 
conical  problem  would  uncover  potential  problem  areas  in  the  computation  of 
these  flows  without  the  added  complexity  of  geometrical  variations. 

The  conical  problem  led  to  the  development  of  special  techniques  to  resolve 
the  vortical  layer  that  inevitably  develops  on  the  body  surface. 

Typical  results  are  shown  for  subsonic  and  supersonic  leading  edge  delta 
wings  and  wing-body  combinations.  These  results  are  also  compared  to  lin¬ 
earized  and  nonlinear  potential  flow  theories  as  well  as  available  experimental 
data. 
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SECTION  I 


INTRODUCTION 

The  standard  assumptions  made  in  order  to  predict  the  flow  about  thin  wings 
at  supersonic  Mach  numbers  and  moderate  angles  of  attack  (i.e.  linearized  po¬ 
tential  flow)  preclude  the  prediction  of  some  of  the  basic  fluid  dynamics  of 
these  flows.  Specifically,  the  formation  of  supercritical  crossflow  regions  in 
the  vicinity  of  the  leeward  wing  leading  edge.  The  assumptions  of  small-per¬ 
turbation  linearized  potential  flow  theory  can  not  predict  supercritical  cross- 
flow  regions,  and  thus  any  embedded  shocks  which  might  result  due  to  recom¬ 
pression.  Consequently,  aerodynamic  design  engineers  who  are  keenly  aware  of 
of  the  shortcomings  and  restrictions  of  applying  linearized  design  and  analysis 
tools  to  supersonic  problems  avoid  possible  nonlinear  areas  by  striving  to  con¬ 
tour  their  wings  at  design  points  (i.e.  lift  coefficients)  that  fall  primarily  with¬ 
in  the  realm  of  linearized  flow  (e.g.  C  ~  .2).  Unfortunately,  this  becomes  an 

JL 

increasingly  difficult  task  as  highly  cambered  wings  are  sometimes  required  to 
achieve  specific  design  and  maneuver  goals.  It  has  been  found  and  noted  in  the 
past  that  optimum  wings  that  have  been  designed  with  linearized  methods  often 
fall  short  of  their  design  goals  (Reference  1)  when  tested  experimentally.  These 
shortcomings  have  often  been  blamed  on  unaccounted-for  viscous  boundary 
layer  effects.  On  the  other  hand,  some  very  important  in  viscid  aspects  of  the 
flow  have  been  neglected  (e.g. ,  crossflow-induced  embedded  shocks).  The  non¬ 
linear  inviscid  behavior  of  the  flow  can  shed  light  on  the  boundary  layer  be¬ 
havior  (i.e.  shock-induced  boundary  layer  separation).  Thus,  with  the  some¬ 
what  unreliable  predictability  of  linearized  methods,  experimental  wind  tunnel 
testing  is  currently  the  final  and  often  expensive  determinant  in  the  design  of 
a  supersonic  wing. 

In  Reference  2  and  3,  a  computer  code  was  developed  primarily  for  comput¬ 
ing  the  three-dimensional  inviscid  hypersonic  flow  over  space  shuttle  configura¬ 
tions  by  integrating  the  nonlinear  Eulerian  equations  of  motion  and  explicitly 


computing  shock  waves.  It  was  shown  that  pressure  and  force  data  could  be 
successfully  predicted  for  this  type  of  vehicle,  thus  significantly  reducing  the 
amount  of  wind  tunnel  testing  necessary  in  the  vehicle  development.  It  is  the 
intent  of  the  present  study  to  apply  procedures  similar  to  those  developed  in 
Reference  3  to  supersonic  vehicles.  The  code  described  in  Reference  3  con¬ 
tained  several  drawbacks  in  that  it  lacked  sufficient  geometric  detail  for  model¬ 
ing  contemporary  complex  supersonic  wing  geometries.  In  addition,  the  pro¬ 
cedure  had  never  been  tested  on  thin-wing  geometries  in  supersonic  flow  where 
typically  very  large  gradients  of  pressure  and  Mach  number  can  be  generated 
in  the  vicinity  of  wing  leading  edges.  An  important  step  in  computing  the  flow 
about  complex  wings  is  to  select  the  appropriate  conformal  mapping  necessary 
to  generate  an  accurate  grid  line  or  computational  mesh  arrangement  in  the 
transversel  plane,  perpendicular  to  the  marching  direction.  This  procedure 
was  first  used  by  Moretti  (Reference  2)  for  three-dimensional  flow  and  was 
shown  to  be  a  significant  factor  in  generating  accurate  and  reliable  results. 
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SECTION  II 


PROBLEM  DEFINITION 

Figure  1  shows  a  sketch  of  the  type  of  flow  field  under  consideration.  Ini¬ 
tially,  conical  wings  will  be  considered  because  of  the  geometrical  simplifica¬ 
tion  of  reducing  the  full  three-dimensional  problem  to  only  two  spatial  coordi¬ 
nates.  In  addition,  the  nature  of  conical  geometries  makes  it  possible  to  study 
the  effect  of  cross  sectional  perturbations  on  the  spanwise  pressure  distribu¬ 
tion.  In  subsequent  sections,  the  conical  solutions  to  the  Euler  equations  will 
be  used  as  starting  conditions  for  three-dimensional  wings. 

The  nature  of  a  conical  flowfield  is  dominated  by  the  “crossflow”  Mach  num¬ 
ber.  The  crossflow  velocity  is  defined  to  be  that  velocity  that  is  tangential  to  a 
spherical  surface  whose  origin  is  the  apex  of  the  conical  wing.  Figure  2  shows 
a  sketch  of  the  cross-sectional  plane  of  a  delta  wing.  At  zero  angle  of  attack, 
the  crossflow  about  a  symmetrical  wing  will  stagnate  at  the  wing  leading  edge 
(see  Figure  2a),  rapidly  accelerate  about  the  upper  and  lower  surfaces  and  then 
stagnate  at  the  upper  and  lower  symmetry  planes.  At  angles  of  attack  greater 
than  zero,  the  crossflow  stagnation  point  moves  slightly  to  the  windward  side  of 
the  leading  edge.  The  resulting  crossflow  undergoes  a  larger  expansion  and 
can  attain  supercritical  values  on  the  upper  surface,  thus  forming  a  supersonic 
crossflow  region  (see  Figure  2b).  The  crossflow  must  stagnate  in  the  leeward 
symmetry  plane  and,  if  the  crossflow  becomes  sufficiently  supercritical,  a 
“crossflow  shock”  is  generated  on  the  upper  surface  or  leeward  side  of  the  wing. 
The  crossflow  shock  facilitates  the  turning  of  the  three-dimensional  velocity 
vector  or  streamline  in  the  direction  of  the  leeward  symmetry  plane,  thus 
satisfying  tho  symmetry  plane  crossflow  stagnation  boundary  condition.  The 
crossflow  shock  is  a  normal  shock  at  the  body  in  order  to  satisfy  the  flow  tan- 
gency  body  boundary  condition,  thus  delineating  the  supercritical  crossflow 
region  from  its  subcritical  counterpart. 


Figure  1  Sketch  of  Supersonic  Fiowfield  About  Elliptic  Cone 
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(«)  ZERO  ANGLE  OF  ATTACK 


SUPERSONIC  CROSSFLOW  (  V  v2p  +  V^,2  >  8) 
CROSSFLOW  SHOCK  ..  ) 


(b)  ANGLES  OF  ATTACK  GREATER  THAN  ZERO 


R  80-1 248-00 1W 


Figure  2  Sketch  of  Crossflow  Stagnation  Streamlines  and  Supersonic  Crossflow  Region 
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The  problem  of  computing  the  external  flow  about  a  supersonic  wing  is 
completely  defined  by  the  wing  geometry,  the  free  stream  Mach  number  and 
angle  of  attack,  and  an  initial  starting  solution  which  satisfies  the  local  boundary 
conditions.  A  finite  difference  marching  technique  can  then  be  applied  to  inte¬ 
grate  the  governing  partial  differential  equations  in  an  assumed  axial  marching 
direction  (Z).  The  only  restriction  is  that  the  axial  Mach  number  must  remain 
supersonic  and  that  the  Courant- Friedrichs- Lewy  (CFL)  stability  condition 
must  be  satisfied  in  marching  from  one  station  to  another. 

In  order  to  apply  the  boundary  condition  at  the  wing  surface  (i.e. ,  vanishing 
of  the  normal  velocity),  the  surface  and  all  its  first  derivatives  must  be  de¬ 
fined.  Therefore  an  analytic  definition  of  the  geometry  is  needed,  with  con¬ 
tinuous  first  derivatives.  The  second  derivatives  of  the  body  geometry  appear 
explicitly  in  the  present  formulation  of  the  problem,  but  continuity  is  not  neces¬ 
sary. 

The  governing  equations  were  defined  in  a  Cartesian  coordinate  system,  thus 
giving  no  special  consideration  to  the  conical  flow  problem.  In  the  marching 
procedure  each  cross  sectional  geometry  is  mapped  into  a  near  circle,  thus 
developing  an  accurate  computational  mesh.  For  the  conical  flow  problem,  the 
mappings  exhibit  conical  similitude  in  the  computational  plane  and  thus  the 
axial  derivatives  will  vanish  if  a  conical  flow  solution  has  been  achieved. 
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SECTION  in 


COMPUTATIONAL  GRID:  CONFORMAL  MAPPINGS 

The  definition  of  a  computational  grid  is  one  of  the  most  crucial  steps  in 
the  building  of  a  numerical  technique.  The  concept  of  using  conformal  mappings 
to  develop  the  computational  mesh  in  three-dimensional  supersonic  flow  prob¬ 
lems  was  originally  proposed  by  Moretti  (Reference  2). 

Three  coordinate  systems  or  spaces  will  be  referred  to:  the  physical  space 
(x,  y,  z),  a  mapped  space  (r,  9,  i  ),  and  a  computational  space  (X,  Y,  Z) ;  see 
Figure  3.  The  coordinate  system  in  the  physical  space  is  Cartesian  and  de¬ 
fines  the  three  velocity  components  that  are  to  be  computed.  The  governing 
equations  are  written  in  the  physical  space  and  then  all  derivatives  are  trans¬ 
formed  into  the  computational  space  where  the  mesh  points  are  at  uniform  in¬ 
tervals  AX,  AY,  AZ. 

Each  j  =  constant  plane  in  the  mapped  space  is  obtained  by  conformally 
mapping  cross  sections  in  the  physical  space  into  near  circles  in  the  mapped 
space.  The  region  bounded  by  the  body  and  the  bow  shock  becomes  an  annulus 
in  the  mapped  space.  The  corresponding  computational  space  is  obtained  by 
normalizing  the  radial  distance  between  body  and  shock  (X-direction)  and  the 
circumferential  distance  between  the  two  symmetry  planes  (Y-direction);  see 
Figure  4.  The  mapped  space  serves  three  purposes.  First,  it  distributes  the 
body  mesh  points  (which  are  evenly  spaced  in  the  computational  plane)  so  that 

\ 

the  necessary  resolution  is  obtained  in  regions  of  large  curvature  where  trunca¬ 
tion  error  may  otherwise  become  too  large.  Second,  the  mapping  makes  the 
body  and  bow  shock  positions  single-value  functions,  r  =  b  (9,  3)  and 
r  =  C(0,a)  in  the  mapped  plane.  Third,  since  shockwaves  embedded  in  the  flow- 
field  become  mesh  lines,  an  accurate  and  stable  computation  is  contingent 
upon  the  shape  of  the  crossflow  shock  not  differing  significantly  from  a  9  = 
constant  mesh  line.  This  enables  one  to  define  a  number  of  regions  in  the 
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computational  space  when  more  than  one  shock  exists  in  a  cross  section.  The 
dotted  lines  in  Figure  4  are  extensions  of  shocks  to  complete  the  boundaries  of 
these  regions.  The  points  on  either  side  of  these  portions  of  the  boundaries 
are  allowed  to  pass  information  across  the  boundary.  The  crossflow  surfaces 
(i.e. ,  extensions  of  crossflow  shocks)  are  defined  as  0  =  constant.  The  con¬ 
stant  taken  is  the  value  of  0  at  the  last  radial  point  on  the  shock. 

Now,  the  transformation  (r,  0,  3)  — (x,  y,  z)  is  considered.  The  mapping 
that  has  been  found  suitable  for  cambered  wings  is  simply  a  single  Karman- 
Trefftz  transformation,  which  may  be  written: 


where  £  =  reie  (mapped  space  variable)  and  W  =  x  +  iy  (physical  space  variable). 
W*  is  the  position  of  the  singularity  of  the  mapping  in  the  physical  spree,  W* 

=  x*  +  iy*  (Figure  5).  The  singularity  of  the  mapping  must  be  placed  near  the 
wing  leading  edge  in  order  to  map  the  wing  cross  section  into  a  near  circle. 

The  singularity  should  be  on  the  wing  cross  sectional  camber  or  mean  line;  this 
condition  defines  y*  and 


x*  =  vXT2  -  RJ;  (2) 

where  xT  =  coordinate  of  the  wing  leading  edge  (Figure  5)  and  Rc  is  the  radius  of 
curvature  of  the  wing  leading  edge.  The  power  6/n  of  equation  (1)  is,  in  gen¬ 
eral,  the  external  angle  of  the  corner  of  which  the  mapping  is  intended  to  eli¬ 


minate  . 


In  the  case  of  a  blunt  wing,  6  =  2n,  and  therefore  equation  (1)  becomes 


In  order  to  compute  W  from  f  (i.e. ,  to  invert  this  mapping)  care  must  be  taken 
to  choose  the  proper  root  of 
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Figure  5  Conformal  Mappings  for  Cambered  Wing 
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rx  =  (uux  + wx)/r  (a) 

r7=  (uu7+  wy)/r  (b) 

6x  =  (uvx-vux)/r2  (c) 

ey=(uvy-vuT)/r2  (d) 


(7) 
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For  the  3  derivatives  of  the  mappings,  the  following  results  are  obtained: 


w  -w*  -(w-w*)  +  wj(w  +  W*) 
»  »  (W*  +  W*)  w*  +  w* 


x4  =Re(Wj) 

y4  =Im(Wa)  (8) 

where 


Wj=xJ+iy* 
w t  =  *l  -iy  4 

Where  xj  =  x*,  y*=y*  since  z=  3  andz4=l.  Now,  since  r,  6,  3  are  inde¬ 
pendent  coordinates,  it  follows  that  r4  =  =  0 .  Therefore: 


ra=r*xi  +  rrya  +  r*zj  =  0 

and 

ei  =  exXi+Ory»  +  ezzi  =  0 


Solving  for  rz  and  6Z,  the  following  results  are  obtained: 

rz  =  -<ryy»  +  rxx4)  (a) 

ez=- (ery»  +9xxS )  (b> 


(9) 


The  second  derivatives  of  this  transformation,  needed  for  the  body  point 
calculation,  are  derived  in  a  similar  fashion. 


The  singularity  of  this  mapping  is  located  inside  the  body  so  that  the  map¬ 
pings  are  never  evaluated  at  singular  points.  It  is  not  necessary  that  the  map¬ 
pings  be  conformal  but  it  has  been  found  that  conformal  mappings  give  the  best 
results  in  terms  of  mesh  point  distribution.  It  is  also  not  necessary  that  r,  6, 

3  be  orthogonal  coordinates  and,  in  general,  they  are  not. 

These  mappings  (Eq.  (1))  use  simple  algebraic  expressions  and  their  coef¬ 
ficients  are  defined  explicitly  so  that  transformation  from  one  space  to  another 
takes  a  minimum  of  time.  Considerable  work  has  been  done  to  develop  conformal 
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mappings  that  can  map  arbitrary  cross  sections  into  circles  or  near  circles 
(Reference  4).  These  generalized  mappings  offer  a  greater  flexibility  than  the 
mappings  used  herein,  but  would  require  a  large  increase  in  computational 
time. 

Some  examples  of  the  mesh  generation  in  the  physical  space  are  shown  in 
Figure  6  for  elliptical  and  cambered  wing  cross  sections.  The  mappings  yield 
the  highest  mesh  density  in  the  vicinity  of  the  leading  edge. 

With  the  mapped  plane  completely  defined,  the  transformation  between  the 
computational  space  and  the  mapped  space  (X,  Y,  Z)  —  (r,  6,  3  )  is  required. 
Consider  a  cross  section  (Z  =  3  =  constant)  with  multiple  shock  waves,  e.g., 
(Figure  7)  bow  shock,  wing  shock  and  tail  shock,  plus  two  crossflow  shocks. 
The  computational  plane  is  divided  into  1C  regions  in  the  Y  direction,  and  LC 
regions  in  the  X  direction.  They  are  ordered  as  in  Figure  4.  The  body  is  de¬ 
scribed  by  r=B(Y,  Z)  and  a  wing-type  shock  as  r  =  C£(Y,  Z)  for  4=1-— LC 
(4=  LC  being  the  Dow  shock).  Similarly,  the  crossflow  shocks  are  described 
as  6  =  Ht  (X,  Z)  for  i  =  1  —  IC  +  1  (i  =  1  is  the  bottom  symmetry  plane  and  i  =  IC  +  1 
is  the  top  symmetry  plane) .  These  surfaces  are  shocks  for  some  range  of  X 
and  Y,  and  arbitrary  (dividing)  surfaces  for  other  values  of  X  and  Y,  as  de¬ 
scribed  previously. 

Now  define  LC  +  1  surfaces  such  that: 

Cj  (Y,  Z)  =  B  (Y,  Z) 

C*(Y,  Z)  =  C  4„  i  ( Y,  Z)  (4=2,  3...  LC  +  1) 

The  transformation  to  the  computational  plane  then  can  be  written  as: 

X=(r-C£)/(C£  +  1-C£)  (a) 

Y  =  (0-H1)/(H1+1-H1)  (b)  (10) 


Z=  i 


(c) 
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Figure  7  Cross  Sectional  Shock  Pattern  in  Mapped  and  Computational  Spaces 
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The  coordinates  X,  Y,  Z  are  not  orthogonal.  The  boundary  C£  in  the  mapped 
plane  becomes  X  =  0  in  region  l  and  Cje+1  becomes  X  =  l.  Similarly,  the  IC 
regions  in  the  Y-direction  in  the  computational  plane  are  bounded  by  Y  =  0  and 
Y  =  1. 


Inverting  this  transformation  yields  the  result: 


r  —  X  (Cf  +  (  —  C£)  +  C£ 

(a) 

6  =  Y(Hltl-H1)  +  H1 

(b) 

(11) 

3  =  Z 

(c) 

Again,  the  derivatives  of  this  transformation,  Xr,  Xe,  Xj,  Yr,  Ye  and  Ya  , 
are  needed. 


First 


rX  “  (Q*  ♦  1  ~ 

(a) 

rY*<CYf  +  1  -  CYj)  X+  CYj 

(b) 

r  z  -  (c  ze  +  ,  ~  CZ^  X  +  CZ£ 

(c) 

ex=(HXf  +  1-HXl)Y+HXl 

(d) 

6y=(Hl+1-H1) 

(e) 

ez  =  (HZt+i-HZl)Y+HZi 

(f) 

O 

II 

X 

(g) 

o 

II 

X 

HO 

(h) 

3Z  =  1 

(i) 

The  Jacobian  of  the  transformation  is  defined  as: 


(12) 


Js  6(r,  6,  3  ) 
d(X,  Y,  Z) 


rX  rY  rZ 

9X  0Y  9Z 
3X  bz 
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and  thus 


X 

Xr  J 

'  8(x.  e,  2,) ' 
L  8(X,  Y,  Z)  _ 

(a) 

•< 

*9 

II 

c_H- 

r d(Y,  e,i)' 
[_  9(X,  Y,  Z) 

(b) 

x»-jl 

"  9(r,  X,  d)l 
_8(X,  Y,Z)  J 

(c) 

< 

<D 

ii 

[  8(r.  Y,  3)1 
L8(Xt  Y.Z)  _ 

(d) 

[  8(r,  9,  X)  1 
L  »(X,  Y,  Z)  . 

|  (e) 

II 

>4 

"  8(r,  0,  Y)  I 
_  9(X,  Y,  Z)  J 

(f) 

After  some  algebraic  manipulations  the  above  derivatives  can  be  written  in 
the  following  form: 


Xr=  1.  / [6  +  XD]  6y+  Dj  CyJ  (a) 

Yr=D,Xr  (b) 

Ytf=l./[A  +  YD2  Ax-D2Hxi]  (c) 

Xg  =  D2  Yg  (d) 


x  -  -  z  +  X<5  yD3  +  Czt+  Cy£  D3I  ^ 

*  [fi+XD46y+CytD4] 

Y,=D3  +  D4X*  (f) 

where 


<5  =Ct  +  t  -  C£ 

(a) 

6V=CY£  +  ,~  CYf 

(b) 

6Z=Cz£4|-CZ£ 

(c) 

(d) 

(14) 
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(15) 


*1+l-HXt 

(e) 

iz  —  Hz 

zi  +  i  zi 

(f) 

(YAx-HXl)/A 

(g) 

(X6  y  +  C  Yj)  / 5 

(h) 

(HZl  +YAZ)/A 

(i) 

(Hjq+YA^/A 

0) 

This  transformation  is  a  modification  of  the  one  used  previously  by  Moretti 
to  solve  numerical  problems  (e.g. ,  see  Reference  5).  Singularities  occur 
when  Ct  =  Cjf  t  j  or  Ht  =  Ht  + 1  and  when  J  =  0.  The  former  occurs  when  two  shocks 
intersect.  This  matter  will  be  discussed  in  the  section  on  “Treatment  of 
Shocks’*.  The  latter  case  occurs  when  the  mesh  lines  X  =  constant  and  Y  =  con¬ 
stant  become  parallel  in  the  physical  plane.  This  can  occur  for  certain  loca¬ 
tions  of  crossflow  shocks.  However  this  problem  can  be  overcome  by  either 
modifying  the  conformal  mappings  so  that  the  cross  section  in  the  mapped 
plane  is  “more  circular,  ”  or  by  using  a  crossflow  shock-type  surface  (which 
acts  like  an  extension  to  a  crossflow  shock)  to  control  the  shape  of  the  mesh 
lines. 


All  shocks  are  defined  in  the  mapped  plane  as  r  =  C(0,9)  and  0  =  h(r,  b  ),  so  that 
C(Y,Z)  =C[0(X8,  Y,Z),d  ],H(X,  Z)  =  h[r(X,  Ys,  Z),  b  ]  (where  X8  and  Y8  are  either  0 
or  1)  and  their  derivatives  C  Y,  C  z,  Hx  and  H  z  must  be  calculated.  The  body 
is  defined  in  the  physical  plane  and  an  iterative  procedure  is  needed  to  de¬ 
scribe  the  body  as  r  =  b  (0,  b  )  from  which  B  (Y,  Z)  can  be  computed.  From  the 
derivatives  be,  bJf  c4,  ha  and  hr,  the  calculation  of  B,  C,  Bz,  and  Cz 
proceeds  as  follows. 

Using  the  notation  of  equations  (11,  12  and  13),  define 


cjsb(e,  ») 


C(SCf.|(fi,}) 
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then 


Cx  =cg  [H(X8,  Z)1+1-H(X8,  Z)t] 


(16) 


HXi  =  hri[C(Y8,  Z)£  + 1  -  C  (Y8,  Z)£] 

Where  again  X8  (the  value  of  X  at  the  shock  C£)  and  Y8  (the  value  of  Y  at  the 
shock  Ht)  are  0  or  1. 


Cz £-CJ£  +  C0£ 


(a) 


^z  =  Y  (Hz  (Xa,  Z)1+1-Hz(Xa,Z)1l+Hz  (X8,  Z)t  (b) 
HZt  =  h  8 1  +  hrt  rZ  (C) 

rz  =  X  [Cz  (Y8,  Z)£  +  1-Cz(Y8,  Z)£]  +  Cz(Y8,  Z)£  (d) 


(17) 


At  the  points  X  =X8  and  Y  =  Y8  at  a  cross  section  Z,  these  equations  result  in 
a  set  of  simultaneous  linear  equations  for  Hz  (X8,  Z)  and  Cz  (Y„,  Z): 


c«£[Y"(h»i*i  h*i  h»i1  »j 

CZ(Y8,  Z)  =  ‘  ^  - r-r- r-r- 

e  1 -ce  [Y8(h  -h  )-h  ] 


(18) 


i  ♦  l 


whereas  Hz  (X8,  Z)  can  be  computed  from  equations  (16).  For  all  other  points, 
equation  (16)  is  used  to  compute  Cz  (Y,  Z)  and  HZ(X,  Z). 

Now  the  computational  plane  and  its  boundaries  are  completely  defined  so 
that,  for  any  mesh  point  (X,  Y,  Z)  in  the  computational  plane,  the  corresponding 
point  (x,  y,  z)  in  the  physical  plane  and  all  the  necessary  transformation  de¬ 
rivatives  can  be  computed. 
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SECTION  IV 


COMPUTATION  OF  REGIONS  OF  CONTINUOUS  FLOW 
In  the  physical  plane  (x,  y,  z)  the  Euler  equations  are: 


wPj  +■  y wz  =  -  ( uPx  +  vPy + y  Uj  +  y vy) 

(a) 

wUj  =  -  (ui^  +  vUy  +  T  Px) 

(b) 

wvz  =  -  (uvx  +  vvy  +  T  Py) 

(O 

(19) 

TPZ  +  wwz  =  -  (uwx  +  vwy) 

(d) 

wSz  =  -  (uSx  +  vSy) 

(e) 

where  T  =  T/T0O,  P  =  ln(p/pJ,  S  =  (S  -  Sj/cVoo  and  all  velocities  are  nondimen- 
sionalized  with  respect  to  Vpx/px  (the  barred  quantities  are  dimensional), 
x  =  x/I,  y=y/T,  z'  =  z/r  (fis  an  arbitrary  length). 

The  equation  of  state  for  an  ideal  gas  becomes: 

ln(T)  =  P(YY1) 

(20) 

The  dependent  variables  are  P,  S  and  the  Cartesian  velocity  components 

u,  v,  w  tFigure  1).  Transforming  all  derivatives  to  the  computational  plane, 

the  following  results  are  obtained. 

fx  =  fxXx  +  fyYx  +  fzZx 

(a) 

fy  =  fX  Xy  +  fy  Yy+fZZy 

(b) 

(21) 

fz=fXXz+fy  Yz  +  fZ  Zz 

(c) 

where  f  is  the  vector  (P,  u,  v,  w,  S)  and 

Xx  =  Xrrx  +  X9ex  +  Xj  b  x 

(a) 

Xjr  =  Xr  ry +X90y  +  X  j  by 

(b) 

(22) 

Xj  -  Xr  rz+Xg0z  +Xfbz 

(c) 
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Similar  expressions  ean  be  written  for  the  Y  and  Z  derivatives.  The  deriva¬ 
tives  of  (X,  Y,  Z)  with  respect  to  (r,  6,  &  )  and  (r,  6,  i  )  with  respect  to  (x,  y,  z) 
have  already  been  discussed  in  the  previous  section. 

Combining  equations  (19'  'nd  (19e),  the  following  form  of  the  Euler  equa¬ 
tions  is  obtained  which  are  used  in  the  present  solution. 

Pz  =  (anPx  +  at2ux+ai3vx  +  auwx  +  biipY  +  bt2UY  +  b13vY  +  b14WY)  (a) 

Uz  =  -(a2iPx  +  a22uX  +  b2lPY  +  b22uY^  tt>) 

(23) 

VZ  =  -  (a31  Px  +  a33vx  +  b31  PY  +  b33vY)  (c) 

wz  =  -(a41Px  +  a42ox  +  a43vx  +  a44Wx  +  b4|  PY  +  b42uY  +  b43vY  +  b44wY)  (d) 

Sz  =  —  (a35Sx  +  b5jSY)  (e) 

where  the  coefficients  appearing  in  equations  (23)  are  defined  in  Appendix  A. 

At  a  data  plane,  Z  =  Z0  =  constant,  all  the  quantities  on  the  right  side  of  equa¬ 
tions  (23)  are  known  and  therefore  the  derivative  fz  can  be  computed  and  used 
to  predict  the  dependent  variables  at  Z  =  Z0  +  AZ. 


The  step  size  AZ  in  the  marching  direction  must  satisfy  the  CFL  (Courant- 
Friedrichs-Lewy)  condition  for  stability  (Reference  6).  If  Ax±  are  the  charac¬ 
teristic  slopes  in  the  X,  Z  plane  and  XY±  are  the  characteristic  slopes  in  the 
Y,  Z  plane,  the  stability  criterion  is  written  as  follows: 


AZX+  -  AX/\x*  AZy*=AY/Ay- 

AZx-  =  AXAx-  azy-=ay/V 


(24) 


Each  of  these  quantities  is  evaluated  for  every  mesh  point  at  the  station 
Z  =  Z0>  and  AZ  is  taken  as  70%  of  the  minimum  of  all  of  these  AZ  values. 


A  modified  MacCormack,  two-level,  predictor-corrector  finite-difference 
scheme  (Reference  7)  is  used  to  integrate  equations  (23).  It  can  be  proved  that 
the  MacCormack  scheme  is  accurate  to  second  order  for  a  linear  system  of 
equations,  so  that  the  truncation  error  is  of  the  form: 
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where  £  is  a  length  in  the  physical  plane.  In  regions  where  dH/dl3  is  large, 
the  mappings  tend  to  assure  A£—  0  so  that  the  truncation  error  remains  small, 
while  keeping  the  total  number  of  grid  points  to  a  minimum. 

Equations  (23)  can  be  written  in  the  following  general  form 

fz  =  [A]fx  +  [B]fY  (25) 

where,  as  previously  defined,  f  is  the  vector  (P,  u,  v,  w,  S)  and  [A]  and  [B]  are 
matrices  of  the  coefficients  of  equations  (23).  With  these  equations,  the 
MacCormack  scheme  proceeds  as  follows. 

Level  one: 

fz  (Z  o) =  (A  1  fx  +  f B]  fY  (all  quantities  evaluated  at  Z  0)  (a) 

f  =  f  (Zq) +  fz  (Z0)  AZ  (f  is  the  predicted  value)  (b) 

Level  two:  (26) 

fz  =  [A  ]  fx  +  [B]  fY  (all  dependent  variables  evaluated  (c) 

with  the  predicted  values  and  all 
independent  variables  are 
evaluated  at  Z  =Zq  +  AZ) 

f(Z0+AZ)  =  (  f +  f  (Zj)  +  fz  AZ)/2  (d) 

The  fx  derivatives  are  taken  one  sided  in  the  positive  X-direction  in  level  one, 
the  fY  derivatives  in  the  positive  Y-direction.  For  level  two  the  direction  of 
these  derivatives  is  reversed. 

This  procedure  defines  all  the  dependent  variables  at  interior  points  of  the 
computational  plane  (1<NN<NC(L)  and  1<MM<MC  (I)),  Figure  7.  The  body 
point  calculation  and  the  shock  point  calculation  will  be  discussed  later.  How¬ 
ever,  note  there  that  all  imbedded  shock  points  have  two  mesh  points  associated 


with  them,  one  for  the  low  pressure  side  and  one  for  the  high  pressure  side, 
both  having  the  same  position  in  the  physical  plane  (Figure  7).  The  low  pres¬ 
sure  side  of  all  shocks  is  computed  following  the  MacCormack  scheme  and 
taking  X  and  Y  derivatives  into  the  low  pressure  region  in  both  levels.  The 
low  pressure  side  of  the  bow  shock  ((NN  =  NC(L),  L  =  LC))  is  defined  by  the 
given  free -stream  conditions. 

The  points  on  the  symmetry  planes  (MM  =  1,  1=1  and  MM  =MC  (I),  I  =  IC)  are 
computed  using  the  same  scheme  and  the  symmetry  conditions  Py=Vy  =  wy 
=  Sy  =  0  and  u  =  0.  The  points  on  the  internal  boundaries  that  are  not  shock 
points  are  also  computed  using  the  MacCormack  scheme.  In  level  one,  the 
points  NN  =  1  and  MM  =  1  are  computed,  taking  the  difference  between  NN  =  1,  NN  =  2 
and  MM  =  1,  MM  =  2  for  the  X  and  Y  derivatives,  respectively.  After  level 

one,  quantities  at  NC(L)  and  MC  (L)  are  updated  (i.e.,  f(NC(L),  M)L  =  f(l,  M)L  +  t 
and  f(N,  MC(I))i  =  f(N,  +  In  level  two  the  points  on  the  other  side  of  the 
surfaces,  NN  =  NC(L)  and  MM  =  MC(I),  are  computed  and  afterward  the  points 
NN  =  1  and  MM  =  1  are  updated. 

A  careful  study  of  the  numerical  results  initially  generated  for  thin,  super¬ 
sonic  wings  indicated  that  the  gradients  in  the  Cartesian  velocity  components  u 
and  v  were  very  large,  while  the  gradient  in  the  velocity  component  in  the  di¬ 
rections  of  the  computational  space  coordinates  X  and  Y  were  much  smaller. 

These  velocity  components  are  defined  by: 

u  =  ulj  +  vly  +  wlj 

v  =  uJx  +  v  Jy  +  wJz 

whereT=VX/|  VXI  and  J*=  VY/I  VY  I 

Solving  for  the  two  Cartesian  velocity  components  u  and  v, 


DET 
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and 


DET 


where 

DET  =  Ij  Jy  —  Jj  Iy 

By  analytically  differentiating  equations  (27)  and  (28),  the  derivatives  of  the 
cartesian  velocity  components  can  be  determined  as  a  function  of  Q,  $  and  their 
derivatives,  which  are  evaluated  numerically.  It  was  found  that  numerically 
differencing  (1  and  $,  instead  of  u  and  v,  reduced  the  truncation  error  of  the 
total  computation  significantly. 

It  should  be  noted  that  Cl  and  $  are  the  normal  and  tangential  crossflow 
velocity  components,  respectively,  on  the  surface  of  the  body.  Thus,  u  is 
constrained  to  vanish  on  the  surface  due  to  the  flow  tangency  boundary  condi¬ 
tion.  A  comparison  of  the  two  sets  of  velocity  components  and  their  variation 
on  the  surface  of  an  ellipse  is  shown  in  Figure  8.  The  improvements  in  the 
computed  results  are  demonstrated  in  Figure  9  at  zero  angle  of  attack.  A 
marked  improvement  in  the  results  at  angle  of  attack  is  demonstrated  in 
Figure  10. 

Special  difference  formulas  must  be  used  to  avoid  difficulties  associated 
with  entropy  discontinuities  (i.e.  ,  vortical  singularities).  These  procedures 
will  be  discussed  in  a  subsequent  section. 

The  calculation  of  interior  points  is  the  most  time  consuming  of  this  com¬ 
putation  mainly  because  it  is  done  many  times.  The  scheme  used  here  keeps 
the  computational  time  to  a  minimum  by  keeping  the  total  number  of  mesh 
points  as  small  as  possible. 
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SECTION  V 


BODY  POINT  COMPUTATION 

The  boundary  condition  at  the  vehicle  surface  is  u  =  0,  where  u  is  the  ve¬ 
locity  normal  to  the  body.  The  entropy  at  the  body  is  computed  by  using  equa¬ 
tion  (23e),  as  for  any  other  mesh  point.  At  the  body,  the  coefficient  of  Sx  (i.e. , 
a 55)  in  equation  (23e)  vanishes  so  that  this  derivative  does  not  affect  the  calcula¬ 
tion  of  entropy  on  the  body. 


To  compute  the  pressure  on  the  body,  the  continuity  and  three  momentum 
equations  (23a,  b,  c,  and  d)  and  the  body  boundary  condition  are  combined  to 
write  a  compatibility  equation  along  the  characteristic  (in  the  X,  Z  plane)  reach¬ 
ing  the  wall  from  the  flow  field. 


where: 


PZ  +  allPX  +  a12uX  +  a13vX  +  a14wX  ~  R1 
uz  +  a21Px  +  a22ux  =  R2 

VZ  +  a31PX  +  a33vX  =  R3 
WZ  +  a41PX  +  a42uX  +  a43VX  +  a44WX  =  R4 


Substituting  for  a12, 


Ri  =  -  (bi,PY  +  b12uY  +  b,3vY  +  b14wY) 
R2  =  —  (b2JPY  "f"  b22uY) 

R3  =  —  ^31RY  +  ^33VY^ 

R4  =  — (b44PY  +  b42  uy  +  b43uY  +  b44wY) 

a13  and  a14  equation  (29a)  becomes 

„  w2yXx  _  .  w2>Xy 

pz  +  anpx  +  a  L  Tx  +  T-2"  ox  _  Ri 


(a) 

(b) 

(29) 

(c) 

(d) 


(30) 


where: 


t  =  u/w  and  <7  =  v/w 


as* 
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Taking  the  difference  between  the  product  of  w  with  equation  (29b)  and  u  with 
equation  (29b)  and  substituting  for  a33,  a43,  a42  and  a44,  the  following  result  is 
obtained. 


(X  TW^  \ 
Xz  +  -  ”  Xyd  j  ' 


(31) 


w^ryT 

+  — —  Xy  rrx  =  wR2  -  uR4 


Similarly,  taking  the  difference  between  the  product  of  w  with  equation  (29c)  and 
v  with  equation  (29d)  and  substituting  for  a33,  a43,  a42  and  a44,  the  following  re¬ 
sult  is  obtained. 

w2crz  +  Px<wa3i  “  va4i>  +  +  XrAj"  +  X*TVx 

X  '  (32) 

w2ctVT  „  _ 

+  ~ - X  ,TX  =  wR3  -  vR4 


The  body  boundary  condition  is  u  =  0.  Since  X  =  constant  is  the  body,  this 
boundary  condition  can  be  written  as 

u  =  uXx  +  vXy  +  wX,  =  0 


or: 


TXj  +  OXy  +  Xj  = 


133) 


Combining  equations  (30),  (31),  and  (32)  and  using  equation  (33),  the  equation 
for  the  characteristic  slopes  can  be  written  in  the  form 


*  w 


2i^lir(YX,  +  u^tx,  ,1 

W  A !  |_  w 1  y  1 


(34) 


where  A_  is  the  slope  of  the  characteristic  reaching  the  wall  from  the  flow  field. 
The  compatibility  condition  along  this  characteristic  is: 


[  A  ”  (UX*  +  vX*>]  <Pz^*px> 


2 

+  7f«,TZ+X,l7z)  + 


^(XlTx+Xy„x)=f 


(35) 
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where: 


H=(X-^<UX-  +  VX'))R' 

+  (wR2  -  uR4)  +  (wR3  -  uR4) 
a1  a1 

The  equation  for  the  body  can  be  written  in  the  form: 

F  =  r  —  B(  Y,  Z) 

Thus  the  body  boundary  condition  is: 

tFx  +  aFy  +  Fz  =  0 

This  equation  holds  for  all  values  of  Z,  thus: 

cTzFy  +  tzFx  =  -(FzZ  +  aFyz  +  tFxZ)  =  R 
Now  using  equation  (10a),  equation  (36)  becomes: 

F  =  r  -  B  =  (C£  -  B)  X 
Thus,  differentiating  equation  (39): 


(36) 


(37) 


(38) 


(39) 


Fx  =  ((C£)x  —  Bx)X  +  (C£  -B)Xx 

Fy  =  [(C£)y  -  By]  X  +  (C£  -  B)  Xy 

and,  for  X  =  0: 

Fx  =  (C£  -  B)  Xx 
Fy  =  (C£  -B)Xy 


(a) 

(40) 

(b) 


(a) 

(41) 

(b) 


Using  equations  (38)  and  (41)  in  equation  (35),  the  following  result  is  obtained: 


R  yw2 

_ A  i  . 


1  ACc  -B)+MX»Tx+X^gx). 

\  ^  (tXx  +  o-Xy) 


- KP x 


This  equation  is  integrated  with  the  same  scheme  used  for  interior  points  with 
the  X-derivatives,  computed  using  three-point  end  differencing  away  from  the 
body. 
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To  compute  the  velocity  components  on  the  body,  an  intrinsic  frame  (I,  J,  K) 
is  used  with  velocity  components  (u,  v,  w).  The  vector  I  is  the  unit  normal  to 

a  a 

the  body,  with  J  and  K  defined  as  follows: 


I  =  Iji  +  I2j  +  13k 

(a) 

J  =  (lx  k)/  |lx  k  |  =  J,i  +  J2j  +  J3k 

(b) 

(43) 

K  =  fx  K  =  Kti  +  K2j  +  K3k 

(c) 

where  i,  j  and  k  are  defined  in  Figure  4.  The  x  and  y  momentum  equations  are 
used  to  compute  v  as  follows.  Equations  (29b)  and  (29c)  are  integrated  using 
the  MacCormack  scheme  to  obtain  u  and  v  and  then  v  is  obtained  from  the 


equation 


v  =  uj,  +  vJ2  (44) 

From  the  integrated  form  of  the  energy  equation,  the  w  component  of  velocity 
can  be  obtained: 


where  T  is  computed  from  P  and  S  and  the  equation  of  state  (19).  The  three  Car¬ 
tesian  velocity  components  are: 


u  =  vJj  +  wKt 

(a) 

v  =  vJ2  +  wK2 

(b) 

(46) 

w  =  vj3  +  wK3 

(c) 
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SECTION  VI 


TREATMENT  OF  SHOCKS 

In  this  section,  the  computation  of  grid  points  on  the  high  pressure  side  of 
all  shock  waves  and  the  detection  of  embedded  shocks  will  be  discussed.  In 
this  area,  we  draw  systematically  from  Moretti’s  extensive  research  on  the 
treatment  of  shock  waves  (Reference  8  and  9). 

The  bow  shock,  all  wing-type  shocks  (i.e. ,  embedded  shocks  which,  in  gen¬ 
eral,  originate  near  the  body  and  move  toward  the  bow  shock  and  are  caused  by 
canopies,  wings  and  vertical  tails,  see  Figure  11)  and  crossflow  shocks  in  the 
flowfield  are  computed  as  discontinuities  satisfying  the  Rankine-Hugoniot  con¬ 
ditions.  The  bow  shock  and  wing-type  shocks  are  defined  in  the  mapped  plane 
by  r  =  c (6,b)  and  the  crossflow  shocks  by  6  =  h(r,a).  At  a  data  plane  Z0,  all 
dependent  variables  are  known  and,  in  addition,  the  quantities  c,  ce,  ca,  h,  hr, 
h  j  are  also  known.  Using  equations  (16  to  18)  C(Y,  Z0),  CY,  Cz,  H(X,  Z0),  Hx 
and  Hz  are  computed  for  all  the  shocks  in  the  flowfield.  At  Z0  +  AZ,  the  posi¬ 
tions  of  the  shock  points  can  be  computed  by  using: 

C(Y,Z0  +  AZ)  =  C(Y,Zl))  +  Cz(Y,Z0)AZ  (a) 

(47) 

H(X,  Z0  +  AZ)=H(X,Z0) +HZ(X,  Z0)AZ  (b) 

Then  Cy(Y,  Z0  +  AZ)  and  Hx(X,Z0  +  AZ)  are  computed  using  central  differences. 
With  these  quantities  and  equations  (16)  c,  ce,  h  and  hr  can  be  computed  at 
Z  =  Z0+AZ. 

After  the  first  level  of  the  MacCormack  scheme,  the  predicted  values  of  the 
dependent  variables  on  the  low  pressure  side  of  all  shocks  are  computed  (the 
variables  on  the  low  pressure  side  of  the  bow  shock  being  the  constant  free 
stream  values).  With  the  low  pressure  side  of  the  shocks  known,  the  high  pres¬ 
sure  side  is  computed  by  an  iterative  process. 


A  value  of  h  4  or  c,  is  guessed,  between  the  values  corresponding  to  an 
infinitely  weak  shock  and  the  value  which  gives  a  subsonic  axial  Mach  number. 
Once  this  guess  is  made,  the  normal  to  the  shock  can  be  defined.  Let 

F=  r  -  c(d ,  3) 


or 

F=  8  — h(r,  3) 

Then  the  normal  to  the  shock  I  is  given  by: 

I  =  (Fzi  +  Fj  +  Fzk)/ V  Fj;  +  F^  +  F|  =  Iji  +  I2j  +  I3k 

where: 

Fx  =  Frrx  +  Fe6x  +  Fa3x 
Fy  =  F^y  +  Fgdy  +  Fj3y 
F2  =  Frrz  +  Fg0z  +  F43  z 

With  the  normal  to  the  shock  defined,  the  Rankine-Hugoniot  conditions  can  be 
applied.  Using  the  subscripts  1  for  the  low  pressure  side  and  2  for  the  high 
pressure  side,  we  have: 

\  =  V,  *  I  (a) 

Mn,  =  Vn/vrf7  (b) 

P2/ Pi  =  (T+  l)/2.  ]/[  1  +Mj^(y  -  l)/2.  ]  (c)  (48) 

P2/P1  =  t(r  +  i)/{r-i)p2/Pi  —  1. ]/[(y+  i)/(y-i)-p2/p,]  (d) 

T2/Tt  =  [T,(p2/p1)]/(p2/pI)  (e) 

vT  =  vT  —  vn  i  (f) 

2  1  1 

where  Mn  and  Vn  are  the  Mach  number  and  velocity  normal  to  the  shock  and  VT 
is  the  velocity  tangent  to  the  shock. 
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An  intrinsic  coordinate  system  is  defined  at  the  shock  with  the  three  direc¬ 
tions  (I,  J,  K),  coordinates  (|,  r),  w)  and  velocities  (u,  v,  w)  such  that  I  is  normal 
to  the  shock  and: 

K  =  (lx  k)/ 1  ix  k  |  =  Kjf  +  kJ  +  K3k 
j  =  lx  K  =  Jti  +  J2j  +  J3k 

In  the  w  plane,  the  characteristic  that  intersects  the  shock  from  the  high 
pressure  side  has  a  slope: 


d|_  _  .  _  (uw  +  a\^Q2+  -a2*) 

dw  ( w2  -  a^) 


(49) 


A  point  (*)  in  the  Z0  data  plane  is  where  the  characteristic  originates.  The 
characteristic  slope  at  the  shock  point  is  first  evaluated  and  then  the  position 
of  the  (*)  point  is  computed  using  the  relations: 


w*  =  -AZ/(K3  +  XI3)  ' 
=  Aw* 

x*  =  XSH  +  +  UJ*K1 

y*  =ySH  +z*h  +  w*k2 


(a) 

(b) 

(c) 

(d) 


(50) 


where  the  subscript  SH  refers  to  quantities  at  the  shock  at  Z0  +  AZ.  Depen¬ 
dent  variables  at  the  (*)  point  are  obtained  by  linear  interpolation.  A  value  of 
the  pressure  on  the  high  pressure  side  of  the  shock  is  computed  using  the  com¬ 
patibility  equation  along  the  characteristic: 


0  =yw2/(aV  u2  +  W2  -  a2  ] 

A  =  (uw  +  aVu2  +  w2  -  a2)/(w2  -  a2) 

R  _  [(u-Aw)  (vP„  +  yvj-  (7  v  uj  yv  w„] 
a/u2  +  w2  -  a2 

dT  =  u*/w*  -ush/wsh 
PgH  =  P*  +  Rw*  -  /3dr 


P2  =  e 


PSH 


(a) 

(b) 

(c) 

(d) 

(e) 

(f) 


(51) 
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where  ft,  X,  and  R  are  averaged  between  the  *  point  and  the  shock  point.  The 
iteration  is  continued  until  this  value  of  pressure  agrees  with  that  computed 
from  the  Rankine-Hugoniot  conditions  for  some  value  of  c  a  or  h  4  .  For  weak 
shocks,  this  iteration  may  converge  to  a  pressure  ratio  P2/p,  <  1;  in  these 
cases  the  value  of  c  a  or  h  a  which  gives  P2/P1  =  1  is  taken. 

Crossflow  shock  points  at  the  body  must  satisfy  the  body  boundary  condition, 
that  is,  the  velocity  normal  to  the  body  on  the  high  pressure  side  of  the  shock 
must  vanish.  This  implies  that  the  shock  normal,  at  the  body,  must  be  perpen¬ 
dicular  to  the  body  normal.  This  condition  gives  a  relationship  between  h  a  and 
hr  at  the  body: 

hr  =  (Fh  —  haFB,)/Fr  (52) 

whe  re : 

FH  -  FBX®X  +  FB]r®r  +  FBz®z 
Fr  =  FBxrx  +  FBrry  +  FBzrz 

and  FBx,  FBy,  FBz  are  x,  y,  z  derivatives  of  FB  =  r  -  b(0,3>  and  r  =  b(0 , 3 )  de¬ 
fines  the  body. 

After  the  second  level  of  the  MacCormack  scheme,  the  corrected,  final  values 
of  the  dependent  variables  on  the  low  pressure  side  of  the  imbedded  shocks,  the 
values  of  Ca  and  ha  computed  after  the  first  level,  and  the  Rankine-Hugoniot 
conditions  are  used  to  compute  the  final  values  of  the  dependent  variables  on  the 
high  pressure  side  of  the  shocks. 

The  first  problem  encountered  when  one  treats  imbedded  shocks  as  discon¬ 
tinuities  is  their  detection.  There  have  been  a  number  of  techniques  proposed 
(see  Reference  9).  One  of  the  earliest  procedures  has  been  found  to  be  well 
suited  for  the  type  of  shocks  encountered  in  this  problem. 

Crossflow  shocks  and  wing-type  shocks  are  detected  in  very  similar  ways. 

For  crossflow  shocks,  the  pressure  distribution  is  monitored  in  the  Y-dlrection 
and,  for  wing-type  shocks,  the  pressure  distribution  Is  monitored  in  the  X-dl- 
rection.  At  a  data  plane  Z  =  Zo,  the  maximum  pressure  gradient  Px  for 
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wing-type  shocks  and  PY  for  crossflow  shocks  is  located.  Then  a  third  order 
polynomial  is  fit  through  the  four  mesh  points  adjacent  to  the  maximum  gra¬ 
dient.  This  polynomial  takes  the  form: 

X  =  a0P3  +  ajP2  +  a2P  +  a3  - 

where 

X  =  X  (for  wing-type  shocks) 

X  =  -  Y  (for  crossflow  shocks) 

and  the  coefficients  a0,  alt  a2  and  a3  are  computed  by  matching  the  curve  fit  to 
the  four  mesh  points.  The  condition  used  to  determine  the  origin  of  a  shock  is 
dx/dP  =  0,  which  implies  dP/dx— •00.  Applying  this  condition  to  equation  (24) 
yields  an  equation  for  Pf  of  the  form: 

d  -  ~~  ai  ^  ~  3ana? 


When  A2  —  3a0a2  =  0,  this  equation  has  one  real  root.  When  this  condition 
is  satisfied,  a  shock  is  inserted  in  the  flow  field  at  Xf  =  a0P3  +  atP^  +  a2Pf  +  a3. 

Crossflow  shocks  are  assumed  to  originate  on  the  body,  so  that  the 
pressure  distribution  is  monitored  in  the  Y -direction  on  X  =  0  (the  body). 

Once  a  shock  is  found  on  the  body,  monitoring  is  begun  at  increasing  values 
of  X  =  constant. 

In  general,  it  is  not  known  at  what  value  of  Y  the  first  shock  point  on  a 
wing-type  shock  will  be  found,  so  that  the  maximum  pressure  gradient  Px  at 
all  values  of  Y  must  be  tested  until  the  first  shock  point  is  detected.  Once  a 
wing-type  shock  is  detected,  additional  shock  points  are  sought  at  the  grid 
points  adjacent  to  the  end  shock  points. 


SECTION  VII 


STARTING  SOLUTIONS 

During  the  course  of  this  study,  several  techniques  were  used  to  “start”  the 
solution  to  the  Euler  equations.  Early  in  the  study,  approximate  Euler  solu¬ 
tions  for  a  circular  cone  at  a  small  angle  of  attack  (References  10  and  11)  were 
employed  to  start  the  computation.  The  geometry  was  then  deformed  continu¬ 
ously  to  the  desired  conical  wing. 

Later  in  the  study,  linearized  potential  flow  solutions  for  delta  wings  with 
subsonic  leading  edges  were  used  and  permitted  the  computation  to  start  with 
the  desired  geometry.  In  the  marching  technique,  the  geometry  was  held  in¬ 
variant.  These  linearized  potential  flow  solutions  are  described  in  the  follow¬ 
ing  section. 

Another  technique  that  was  employed,  for  both  conical  and  three  dimensional 
wings,  was  to  use  an  exact  Euler  solution  to  start  the  computation.  For  these 
cases,  the  geometry  would  also  be  continuously  deformed.  For  example,  a 
conical  Euler  solution  would  be  used  to  start  the  computation  for  another  coni¬ 
cal  wing.  The  geometry  would  start  conically  then  be  continuously  deformed  to 
another  conical  geometry. 

1.  LINEARIZED  POTENTIAL  FLOW  THEORY 

Since  the  supersonic  flowfield  about  thin  wings  is  of  primary  interest,  a 
starting  solution  more  representative  of  this  type  of  flowfield  was  sought.  As  a 
first  step  in  this  direction,  the  linearized  flow  solution  for  the  flow  about  thin 
elliptic  cones  with  subsonic  leading  edge  at  zero  angle  of  attack  was  used.  This 
solution  can  be  easily  found  in  Reference  8  for  linearized  boundary  conditions 
applied  on  the  plane  Z  =  0. 

The  general  solution  for  the  entire  Cartesian  velocity  field  about  the  elliptic 
cone  can  be  written  in  complex  form  as 
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u  =  Re 


-  nr 


1  [  r  tan©  +  F(0  , K)  -E(0,K)] 

1/2  \ 


m2  f  1  —  t  2  1; 
1  -  m2  1_  T  5  -  m2J 


|1  -  m‘ 

f 

v  =  Re 

V" 

w  =  Re  |  -  m2/(  m2  -  t2]  1/2 1 
where  F  and  E  are  incomplete  elliptic  integrals  of  the  first  and  second  kind. 
Also 

K  =V  1  -  it? 


(53) 
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.  -i  ri-T2 11/5 
=sin  [t^\ 


where  t  is  a  complex  variable  defined  by  the  transformations 

2e 


T  = 


l+€^ 


€  = 


y  +  iz 


X  +  Vx2  —  y2  —  z2 
and  m  denotes  the  leading  edge  value  of  t. 


With  these  transformations,  the  entire  flowfield  about  the  elliptic  cone  can 
be  defined  up  to  the  Mach  cone  where  all  disturbances  vanish.  Special  considera¬ 
tion  must  be  given  to  these  expressions  for  the  elliptic  integrals  in  the  symmetry 
planes  Z  =  0  and  y  =  0  because  one  or  the  other  of  the  incomplete  elliptic  inte¬ 
grals  reduce  to  the  complete  integral, 

A  computer  subroutine  was  coded  to  compute  this  flowfield  to  be  used  as  a 
starting  solution  for  the  nonlinear  code.  Several  problems  arose  in  trying  to 
implement  the  solutions  because  of  the  singularities  that  exist  at  the  leading 
edge  due  to  the  linearized  boundary  conditions.  As  a  result,  leading  edge  cor¬ 
rections  had  to  be  fitted  to  the  linearized  solution  to  obtain  finite  values  of  ve¬ 
locity  and  pressure  at  the  leading  edge  (see  Figure  12).  The  resulting  solution 
then  had  to  be  modified  in  order  to  satisfy  the  exact  nonlinear  boundary  condi¬ 
tions.  The  exact  equation  for  the  pressure  also  had  to  be  used  to  be  consistent 
with  the  nonlinear  equations.  The  subroutine  was  then  made  compatible  with  the 
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Figure  12  Elliptical  Cone  Thickness  Solution,  Subsonic  Leading  Edge 


nonlinear  code  for  its  starting  solution.  The  shock  was  initially  defined  as 
the  free-stream  Mach  cone  about  the  elliptic  cone. 


In  conjunction  with  the  thickness  solution,  a  lifting  solution  for  a  flat  plate 
delta  wing  with  subsonic  leading  edges  can  be  found  in  Reference  12.  The  gen¬ 
eral  solution  for  the  entire  Cartesian  velocity  field  from  the  body  to  the  Mach 
cone  can  also  be  expressed  in  complex  form  as, 


This  solution  in  linear  combination  with  the  solution  for  the  velocity  field 
about  elliptic  cones  was  used  to  represent  an  initial  starting  solution  for  the 
flow  about  a  lifting  delta  wing.  The  Mach  cone  rotated  at  the  angle  of  attack  was 
used  as  an  initial  guess  for  the  location  and  shape  of  the  bow  shock.  In  addi¬ 
tion,  leading  edge  corrections  also  were  implemented  to  this  solution  to  eliminate 
problems  with  the  leading  edge  singularity  associated  with  the  linearized  poten¬ 
tial  flow  solution.  To  be  consistent  with  Euler’s  equations,  the  exact  expres¬ 
sion  for  the  pressure  coefficient  was  used  and  the  entire  initial  incremental  en¬ 
tropy  field  was  set  to  zero. 

The  useful  aspect  of  employing  a  linearized  starting  solution  is  that  lower 
free-stream  Mach  numbers  and  larger  delta  angles  can  be  achieved  which  were 
not  possible  using  the  cone  as  a  starting  solution.  In  addition,  as  mentioned 
earlier,  the  entire  entropy  field  developed  gradually  from  the  irrotational 
starting  solution.  This  allowed  for  a  better  understanding  of  the  development 
of  the  vortical  singularities  and  their  effect  on  the  solution. 
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SECTION  VIII 

SPECIAL  TECHNIQUES  FOR  CONICAL  FLOW 
1.  CROSSFLOW  SHOCK  SHAPE 

Several  problems  arose  in  the  detection  of  crossflow  shocks  and  in  the 
transient  shape  of  these  shocks  during  the  iterative  marching  technique  used 
to  obtain  a  conical  solution  to  Euler’s  equations. 

Figure  13  shows  a  schematic  of  a  typical  supersonic  crossflow  region. 

Since  a  crossflow  shock  can  be  detected  anywhere  along  the  sonic  line  between 
the  symmetry  plane  and  supersonic  crossflow  region,  a  shape  criterion  for  the 
crossflow  shock  had  to  be  imposed  to  ensure  a  stable  computation. 

The  crossflow  shock  is  located  along  the  sonic  line  and  its  strength  decays 
radially  along  this  line.  As  the  shock  becomes  weak,  its  shape  basically  follows 
the  shape  of  the  sonic  line.  The  computational  scheme  assumes  the  shock  is  an 
M-grid  line  boundary  (see  Figure  14a).  If  weak  crossflow  shock  points  are  de¬ 
tected  whose  shape  significantly  differs  from  the  shape  of  an  original  M-grid 
line  developed  from  the  mapping,  the  original  near  orthogonality  of  the  two- 
dimensional  mesh  will  be  disturbed,  leading  to  an  unstable  computation  and 
spurious  numerical  results.  Thus,  a  shape  criterion  was  imposed  to  eliminate 
this  problem  and  retain  a  reasonable  grid  line  arrangement.  Crossflow  shock 
shapes  as  shown  in  Fig.  14b  promote  a  more  stable  computation. 

During  the  iterative  marching  to  a  conical  solution,  certain  transient  and 
sometimes  unstable  phenomena  occurred  which  were  all  associated  with  the 
vortical  layer  that  builds  upon  the  body  surface.  It  was  observed  that  the  cross- 
flow  shock  point  on  the  body  would  become  displaced  relative  to  the  field  shock 
points,  as  demonstrated  by  Figure  15a.  Normally,  the  radial  shock  derivative 
at  all  field  shock  points  is  determined  numerically  using  central  difference 
formulae.  It  was  found  that  this  procedure,  when  implemented  early  during 
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Figure  13  Grid  Line  Configuration  Relative  to  Supercritical  Region 
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Figure  14  Detection  of  Stable  Crossflow  Shock  Shape 


(a)  TRANSIENT  SHOCK  PRIOR  TO  CONVERGENCE  I  ,b)  INSTABI LITV  OF  SHOCK  DURING 
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Figure  15  Shock  Shapes  During  Iterative  Process 
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the  convergence  process,  could  lead  to  an  instability  in  the  shoek  shape1  as 
demonstrated  by  Figure  15b.  The  field  shoek  position  would  beeome  unstable 
and  travel  into  the  supersonie  region,  beeome  weak,  and  eventually  fail  to  con¬ 
verge.  It  was  determined  later  in  the  study  that  this  problem  was  most  prob¬ 
ably  eaused  by  the  buildup  of  a  vortieal  layer  on  the  body  surfaee.  This  vorti- 
eal  layer  eannot  be  resolved  numerieally  without  speeial  consideration.  Thus, 
it  was  necessary  during  the  transient,  while  the  shoek  points  were  moving  to  a 
converged  location,  to  detaeh  the  field  shoek  point  derivatives  from  the  location 
of  the  body  shoek  point.  An  assumption  on  the  direction  of  the  shock  normal  in 
the  eross  seetional  plane  had  to  be  prescribed.  The  shoek  normal  in  this  plane 
was  assumed  to  be  tangential  to  the  loeal  N-grid  line  direction. 

2.  VORTICAL  LAYER 

Prior  to  eonvergenee,  eertain  problems  arose  in  the  pressure  and  erossflow 
Maeh  number  distribution  between  the  erossflow  shoek  and  symmetry  plane. 
Severe  oscillations  developed  in  both  pressure  and  crossflow  Maeh  number  in 
this  region  as  soon  as  the  erossflow  shoek  beeame  strong  and  approaehed  eon¬ 
vergenee.  To  examine  the  problem,  a  numerical  seheme  was  devised  to  allow 
the  insertion  of  more  grid  points  in  this  region.  Two  M-grid  regions  were  as¬ 
sumed:  one  being  the  region  between  the  windward  symmetry  plane  and  the 
erossflow  shoek,  and  the  other  being  the  region  between  the  erossflow  shoek  and 
leeward  symmetry  plane.  Separating  the  grid  into  two  regions  allowed  for  dif¬ 
ferent  grid  spaeing  in  the  two  regions  in  the  computational  spaee.  Thus,  to  in¬ 
vestigate  the  problem  between  the  erossflow  shock  and  leeward  symmetry  plane, 
more  points  were  inserted  in  this  region.  Figure  16  shows  a  typieal  distribution 
of  pressure,  erossflow  Maeh  number,  and  entropy  that  would  result  just  prior 
to  eonvergenee.  Inserting  more  points  in  this  region  permitted  resolution  of 
any  steep  gradients.  A  large  gradient  in  entropy  was  revealed.  Near  the  shoek, 
the  entropy  corresponded  to  the  entropy  being  developed  by  the  shoek,  and  at 
the  symmetry  plane  the  entropy  was  still  negligible.  Thus,  a  wave  in  entropy 
was  slowly  traveling  towards  the  leeward  symmetry  plane.  The  steep  gradient 
in  entropy  was  also  eausing  the  oscillations  in  both  pressure  and  erossflow  Maeh 
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number  in  this  region.  Since  the  body  surface  is  a  streamline,  the  entropy  at 
the  shock  must  propagate  to  the  symmetry  plane,  thus  constituting  the  vortical 
singularity.  The  computation  was  allowed  to  proceed  in  the  hope  that  the  en¬ 
tropy  would  propagate  to  the  symmetry  plane,  but  the  process  turned  out  to  be 
too  slow  in  converging  since  the  entropy  propagates  at  the  crossflow  velocity  and 
thus  becomes  an  asymptotic  process.  Thus,  a  new  computational  technique 
would  be  required  to  speed  up  this  process,  allowing  the  vortical  singularity 
to  develop.  In  addition,  a  new  computational  technique  would  be  required  to 
handle  the  vortical  layer,  and  the  discontinuity  and  large  gradients  in  entropy 
that  would  develop  between  the  body  surface  and  the  first  N-grid  line  off  the 
body. 

Figure  17  shows  a  basic  sketch  of  the  vortical  singularities  and  the  entropy 
field  occurring  with  a  crossflow  shock.  The  entropy  on  the  surface  of  the  body 
between  the  windward  symmetry  plane  and  low  pressure  side  of  the  crossflow 
shock  comes  from  the  stagnation  streamline  which  originates  at  the  bow  shock. 
The  entropy  on  the  body  in  Region  II  comes  from  the  crossflow  shock.  A  strong 
vortical  singularity  is  developed  in  the  leeward  symmetry  plane  on  the  body. 
Here,  the  body  surface  is  a  streamline  and  the  entropy  must  be  equal  to  that 
on  the  high  pressure  side  of  the  crossflow'  shock.  At  the  first  field  point  in  the 
leeward  symmetry  plane,  the  entropy  must  be  that  corresponding  to  the  entropy 
of  the  streamline  coming  from  the  bow  shock.  At  a  high  enough  angle  of  attack, 
a  discrete  bow'  shock  does  not  exist  in  the  leeward  symmetry  plane,  and  thus 
the  entropy  is  zero.  Hence,  the  vortical  singularity  corresponds  to  the  jump 
in  entropy  on  the  body  developed  by  the  crossflow  shock  and  the  entropy,  if  any, 
developed  across  the  leeward  plane  bow  shock  point. 

To  resolve  this  singularity  and  to  speed  up  the  process  of  the  entropy  de¬ 
velopment  on  the  leeward  symmetry  plane  from  the  crossflow  shock,  a  new  com¬ 
putational  technique  was  developed.  This  scheme  is  only  imposed  when  the 
crossflow  shock  attains  maximum  strength  and  is  close  to  convergence.  The 
entropy  on  the  body  between  the  leeward  symmetry  plane  and  the  crossflow 
shock  is  set  equal  to  the  entropy  at  the  crossflow  shock.  At  the  first  N  ring 
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off  the  body,  outward  difference  numerical  derivatives  in  the  radial  direction 
are  used.  The  entropy  on  the  body  is  undisturbed  by  the  first  field  ring  since 
the  coefficients  of  its  normal  derivatives  are  zero.  The  entropy  at  the  field 
ring  is  unaffected  by  the  body  entropy  because  of  the  windward  differencing  in 
the  radial  direction  at  this  ring.  All  other  field  points  are  still  computed  using 
central  difference  formulas. 

To  eliminate  the  oscillations  in  pressure  and  cross  flow  developed  as  a  re¬ 
sult  of  the  steep  entropy  gradient  between  the  shock  and  symmetry  plane,  a  one- 
step  numerical  damping  or  artificial  viscosity  was  imposed  in  this  region. 

Some  typical  results  are  shown  in  Figure  18  for  the  leeward  side  of  the  body. 

As  indicated  by  the  figure,  when  this  procedure  was  adopted  the  results  im¬ 
proved  markedly.  Both  the  pressure  and  crossflow  are  smooth  and  show  no 
oscillations.  The  entropy  distribution  shows  a  discrete  jump,  representing  the 
vortical  singularity,  from  the  body  surface  to  the  first  field  ring.  In  addition, 
a  discrete  jump  in  radial  velocity  also  occurs  at  the  vortical  singularity. 


47 


124S-013W 


N  GRID 


Figure  18  Improved  Numerical  Results 
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SECTION  IX 


GEOMETRY 

1.  ANALYTICAL  SPANWISE  CONICAL  DEFINITION 

A  geometry  definition  based  on  simple  conical  parameters  was  developed 
and  was  found  to  be  very  useful  in  defining  both  conical  and  three-dimensional 
wings.  The  cross  section  was  simply  defined  as  the  superposition  of  a  span- 
wise  thickness  and  camber  distribution.  The  cross  section  was  restricted  to 
spanwise  elliptical  thickness  and  circular  arc  camber.  In  a  subsequent  section, 
a  more  flexible  geometry  package  is  described  for  general  three-dimensional 
wings  with  chordwise  definition. 

In  this  analytical  representation  of  the  wing  geometry,  the  geometric  pa¬ 
rameters  are  specified  as  parametric  angles.  Thus,  a  wing  cross  section  can 
be  entirely  specified  by  five  angles:  6F,  6E,  6T,  6TC,  6P.  The  symbols  6F, 
and  6t  represent  planform  angles  defining  the  body  lines  of  the  delta  wing,  6P 
an  angle  controlling  the  centerline  thickness  variation,  and  6TC  the  leading  edge 
camber  angle  in  the  cross  sectional  plane  (see  Figure  19).  Cubic  equations 
are  used  to  curve  fit  the  angles  in  the  general  form: 

6(z)  =  a(z  -  z0)3  +  b(z  -  z0)2  +  c(z  -  z0)  +d 
The  spatial  location  of  the  body  lines  can  then  be  simply  defined  as: 

x  =  z  tan(S(z)) 

The  thickness  of  the  wing  was  assumed  elliptical  but  was  allowed  to  have  an 
inboard  flat  section.  The  camber  line,  defined  separately,  was  also  allowed 
to  have  a  flat  inboard  section  with  a  circular  arc  camber  line  to  the  leading 
edge.  To  define  the  maximum  camber  and  twist  of  the  cross  section,  a  leading 
edge  camber  angle  is  specified  (6TC).  The  thickness  and  camber  flats  are  de¬ 
fined  separately  and,  thus,  do  not  have  to  coincide.  The  wing  cross  sectional 
geometry  was  then  simply  defined  by  adding  and  subtracting  the  thickness  to 
the  camber  line. 
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All  first  derivatives  of  the  surface  could  then  be  defined  analytically  through 
the  surface  normal  vector.  The  entire  conical  wing  was  then  defined  by  three 
continuous  wing  sections.  The  geometry  definition  was  flexible  enough  to  allow 
the  initial  section  to  begin  as  a  circular  cone.  The  cone  was  then  deformed  con¬ 
tinuously  through  two  wing  sections.  The  third  section  represented  the  in¬ 
variant  conical  portion  of  the  wing.  The  parametric  definitions  remained  un¬ 
changed  for  all  three  wing  sections,  with  only  the  coefficients  of  the  cubic 
angle  equations  varying  from  one  section  to  another. 

As  an  example,  for  a  simple  three-dimensional  elliptical  cross  section  wing, 
the  surface  normals  would  be  defined  from  the  parametric  cross  section  equa¬ 
tion, 

y  =  Pt(z)[x|(z)-x2],/2 

where  Pt(z)  is  a  cubic  centerline  thickness  function  and  XT  =  z  tan6T(z)  where 
(5t  is  a  cubic  planform  angle  function  in  Z.  The  surface  in  three  dimensions 
can  then  be  defined  simply  as, 

F (x,  y,  z)  =  y  -f  (x,  z) 
and  the  surface  normal  vector  as, 

n*=  (iy  “  f*  i*  "  fz  i*)/d  +  f  x  +  f  z> 1/2 

Figure  20a  shows  an  example  of  the  body  line  variation  to  achieve  a  cambered 
wing  cross  section.  Figure  20b  shows  an  example  of  the  cross  sectional  de¬ 
formation  of  a  circular  cone  to  an  elliptical  cross  section.  Figure  20c  further 
shows  the  deformation  from  an  elliptical  to  a  cambered  cross  section. 

It  was  a  trivial  task  later  in  the  study  to  add  additional  sections  to  the  wing. 
The  present  capability  of  the  geometry  allows  for  variation  from  a  circular 
cone  to  a  conical  wing,  then  an  additional  variation  to  a  second  conical  wing. 


A  wing-body  capability  was  then  added  to  the  second  conical  wing  section. 
This  was  achieved  by  allowing  for  a  cubic  body  section  that  follows  a  body  line. 
Thus,  the  last  section  of  the  wing  can  be  a  conical  wing  or  wing-body.  Figure 
21  shows  a  sketch  of  the  parameters  that  control  the  body  development.  Thus, 
this  geometry  package  can  be  used  for  conical  wings,  conical  wing -bodies, 
and,  if  done  properly,  for  an  entire  three-dimensional  wing-body  vehicle. 
Figure  22  shows  an  example  of  the  usage  of  this  capability  for  computing  a 
conical  wing-body  solution  commencing  with  a  conical  wing  cross  section. 

2.  SURFACE  PATCH  GEOMETRY 


The  analytical  geometry  programs  described  in  the  previous  section  were 
an  advantageous  device  for  computing  results  for  a  variety  of  wing  shapes 
where  the  geometrical  parameters  could  be  computed  exactly  and  without  nu¬ 
merical  errors.  Unfortunately,  not  all  wing  shapes  fall  into  the  category  of 
these  programs,  and  a  more  flexible  geometry  program  was  desirable. 


A  geometry  program  was  desired  which  could  analytically  simulate  very 
complex  wing  shapes  with  a  minimum  of  errors.  After  some  investigation, 
the  bi-cubic  surface  patch  formulation,  developed  in  Reference  13  and  put  to 
practical  use  in  describing  aerodynamic  vehicles  in  Reference  10,  was  found 
to  be  an  adequate  procedure.  The  equation  for  the  Cartesian  coordinates  of  a 
surface  bi-cubic  patch  can  be  written,  in  terms  of  parametric  variables  u  &  w 
(using  matrix  notation)  as  (see  Figure  23): 


x=  [u3  u2  u  1]  [SJ 
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w3 

w2 

w 

1  . 


(b)  (55) 


z=[u3  u2  u  1)  [Sj) 


(c) 
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Figure  21  Wing-Body  Definition 


where 


Sj  =  M  Bj  M*,  i  =  x,  y,  z 

or,  in  vector  notation, 

V  =  [xy  z]  =  USW* 

The  boundary  condition  matrix  B  is  a  4x  4  matrix  containing  the  four  corner 
point  spatial  and  derivative  information.  The  matrix  M  is  a  constant  matrix 
derived  from  the  desired  blending  functions.  The  computer  code  described  in 
Reference  10  uses  the  patch  theory  developed  in  Reference  12  for  the  geometry 
of  complex  aerodynamic  vehicles.  The  code  uses  the  Harris  geometry  input 
data  format  that  is  commonly  used  by  aerodynamicists.  The  wing  sections  are 
defined  chordwise  for  each  span  location.  Spline  fitting  provides  the  necessary 
information  for  developing  the  boundary  matrices  B  for  each  patch. 

The  computer  code  described  in  Reference  10  was  extensively  modified  and 
incorporated  into  the  finite  difference  code.  The  finite  difference  procedure  re¬ 
quires  the  surface  coordinate  y  and  surface  normal  vector  components  at  a  given 
axial  location  z  and  span  location  x  of  the  wing.  Thus,  with  two  of  the  Cartesian 
coordinates  specified,  equations  (55a)  and  (55c)  yield  two  simultaneous  bi¬ 
cubic  equations  for  u  and  w  that  must  be  solved  iteratively.  The  solution  is 
accomplished  quickly  by  a  Newton  iterative  method  for  two  simultaneous  equa¬ 
tions  and  an  error  extrapolation  scheme. 

Thus,  if 

F  (x,  u,  w)  =  x  -f  (u,  w) 

G(u,  w)  s  z  -  g  (u,  w) 

and  (56) 

U|  ♦  i  =  ut  +  6 

WlM  =W,+€ 
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where 


/  Fi  Fwl\  /Fut  fA 

_  \Gt  Gwl  /  €  _  \G  U|  Gt/ 

8(F,G)t  5  9(F,  G)t 

8(u,w)  8(u,  w) 


and  ult  wt  are  an  initial  guess  or  subsequent  iterations  until  a  convergent 
solution  has  been  attained. 


Determining  u  and  w  explicitly  defines  the  surface  coordinates  y  from  equa¬ 
tion  (55).  The  surface  normals  at  the  point  are  then  determined  by  the  cross 
product  of  the  two  tangent  vectors  given  by 


and 


where 
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In  terms  of  practical  usage,  Figure  24  shows  an  example  of  a  delta  wing 
modeled  using  bi-cublc  surface  patches.  Each  patch  equation  requires  a 
(4x4x3)  matrix  of  boundary  condition  information.  The  core  requirements 
necessary  to  store  all  the  patch  matrices  are  beyond  practicality.  Thus,  the 
patch  matrices  are  written  on  an  external  device  and  only  local  blocks  of  patch 
matrices  are  stored  in  the  vicinity  of  the  axial  station  being  currently  computed. 


324B-01  7W 


Figure  24  Panel  Model  of  Delta  Wing 
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SECTION  X 


DISCUSSION  OF  RESULTS 


1.  CONVERGENCE 

The  axial  variation  of  surface  pressure,  in  the  symmetry  plane  and  at  the 
wing  tip,  is  shown  in  Figure  25;  the  geometry  is  a  6:1  elliptic  cone,  and  the 
free  stream  conditions  are  M00  =  1.97,  and  a  =  0.  The  iteration  was  started  at 
Z  =  0. 1  with  a  circular  cone  solution  of  half  angle  18.39°.  The  geometrical  de¬ 
formation  expands  and  recompresses  the  flowfield  between  Z  =  0.1  and  Z  =  1.0 
and  the  geometry  becomes  conical  at  Z  =  1.2.  The  symmetry  plane  and  tip  pres¬ 
sures  are  both  converged  at  Z  =  4.  Figure  25  shows  that  these  values  continue  to 
remain  invariant.  At  zero  angle  of  attack,  the  two  pressure  levels  are  a  good 
indication  of  the  convergence  of  the  flowfield  on  the  body.  For  these  results,  a 
20  x  25  radial  by  circumferential  grid  was  used.  Figure  26  is  a  comparison  of 
the  computed  spanwise  pressure  distribution  with  the  experimental  data  of  Jor¬ 
gensen,  Reference  7.  In  addition,  the  second  order  result  of  Van  Dyke,  Refer¬ 
ence  15,  is  also  shown. 

At  angles  of  attack  high  enough  to  produce  a  crossflow  shock,  the  shock 
pressure  jump  and  position  were  studied  to  determine  convergence.  Figure 
27a  shows  the  axial  distribution  of  the  surface  pressure  in  the  windward  and 
leeward  symmetry  planes  and  the  high  and  low  pressure  side  of  the  crossflow 
shock.  The  configuration  was  the  same  6:1  ellipse  and  the  free  stream  condi¬ 
tions  were  M.  =  1.97  and  a  =  10°.  Figure  27a  shows  all  these  pressure  levels  to 
be  converged  after  Z  =  3.0.  Figure  27b  shows  the  cross  sectional  surface  pres¬ 
sure  distribution  and  geometry.  For  computations  at  angle  of  attack, 
a  25 x  35  grid  was  generally  used.  The  results  of  the  computations  at  the  higher 
angles  of  attack  (a  >  10°)  exhibited  pressure  oscillations  on  the  high-pressure 
side  of  the  crossflow  shock,  as  shown  in  Figure  27b  and  later  in  this  section. 
These  oscillations  seem  to  be  generated  at  the  vortical  singularity  and  move 
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Figure  25  Axial  Variation  of  Surface  Pressure,  Elliptic  Cone  (a/b  «  6),  =  1.97,  a  =  0 
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Figure  26  Comparison  with  Experimental  Data  and  Second  Order  Theory,  Elliptic  Cone 
(a/b  ■  6),  Mgo  “  1.97,  a -0. 
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toward  the  crossflow  shock  as  the  computation  is  continued.  The  amount  of 
difficulty  encountered  in  the  computation  seems  to  be  directly  related  to  the 
strength  of  the  crossflow  shock.  The  resolution  of  this  difficulty  will  be  demon¬ 
strated  in  subsequent  results. 

Figure  28  shows  the  plan  view  of  the  wing  and  the  position  of  the  cross- 
flow  shock  on  the  body  which  was  first  detected  at  Z  =  0.75.  The  shock  position 
approaches  a  conical  ray  at  Z  =  2.3  and  thereafter  remains  on  that  ray. 

2.  EFFECT  OF  SPANWISE  CAMBER 

As  a  preliminary  investigation  of  the  possible  effects  of  conical  camber, 
the  flowfields  about  several  wings  with  various  forms  of  circular  arc  spanwise 
camber  were  computed.  In  Figure  29,  the  computed  spanwise  pressure  distribu¬ 
tions  are  presented  for  a  symmetrical  elliptical  wing,  a  wing  with  half  semi¬ 
span  circular  arc  camber,  and  a  wing  with  circular  arc  camber  over  its  entire 
semi-span  at  zero  angle  of  attack  and  M«  =  2.0.  Both  camber  lines  terminate 
with  a  20°  leading  edge  camber  angle.  The  conical  camber  effectively  results  in 
leading  edge  droop  and  thus  reduces  the  effective  angle  of  attack  of  the  wing, 
resulting  in  a  corresponding  loss  in  lift.  The  pressure  distribution  of  Figure  29 
confirms  this  expected  result  with  higher  pressures  being  computed  on  the  upper 
surface.  The  pressure  distributions  also  indicate  a  slight  movement  of  the 
stagnation  point  to  the  upper  surface. 

Figure  30  shows  the  computed  pressure  distributions  for  a  wing  with  half 
semi-span  circular  arc  camber  and  leading  edge  camber  angles  of  20°  and  30°. 
Figure  31  shows  the  corresponding  crossflow  Mach  number  distributions.  At 
20°  camber,  a  small  region  of  supercritical  flow  is  developed  on  the  lower  sur¬ 
face  of  the  wing  in  the  vicinity  of  the  leading  edge.  At  30°  camber,  this  region 
becomes  more  extensive,  reaching  a  crossflow  Mach  number  of  1.4  and  causing 
a  crossflow  shock  to  form  on  the  lower  surface  at  a  span  location  10%  inboard 
from  the  wing  leading  edge.  This  type  of  behavior  brought  on  by  rather  large 
camber  angles  can  cause  the  boundary  layer  to  separate  on  the  lower  surface 
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Figure  29  Effect  of  Spanwise  Camber  on  Pressure  Distribution,  M, 


of  the  wing  at  small  angles  of  attack.  Leading  edge  separation  on  the  lower 
surface  was  detected  experimentally  in  Reference  15. 

Figures  32  and  33  show  the  effect  of  full  semi-span  circular  arc  camber 
on  the  spanwise  pressure  distribution  for  two  wings  with  different  cross-sec¬ 
tional  thicknesses.  Figure  32  shows  a  14:1  elliptical  wing  at  a  =  3°.  The  sym¬ 
metrical  wing  develops  a  supercritical  region  on  the  upper  surface  and  forms 
a  small  crossflow  shock  at  the  end  of  a  rather  steep  compression.  The  effect  of 
cambering  the  wing  is  to  completely  eliminate  the  suction  peak  and  correspond¬ 
ing  supercritical  region  in  the  vicinity  of  the  leading  edge.  Figure  33  shows  a 
similar  result  for  an  18:1  elliptical  wing  at  a  =  2. 5°.  The  symmetrical  wing  in 
this  case  develops  a  somewhat  stronger  shock,  which  is  located  closer  to  the 
leading  edge.  Cambering  this  wing  also  completely  eliminates  the  supercritical 
region  near  the  leading  edge.  In  both  cases,  a  loss  in  lift  accompanies  the 
camber  effect.  An  improvement  in  the  drag  characteristics  might  be  expected  as 
a  result  of  a  thrust  acting  on  the  forward-facing  area  of  the  wing,  produced  by 
the  leading  edge  droop.  In  any  event,  an  apparent  delay  in  the  formation  of  the 
crossflow  shock  is  apparent. 

3.  COMPARISON  WITH  LINEARIZED  POTENTIAL  FLOW  THEORY 

The  failure  of  linearized  theory  in  the  vicinity  of  the  leading  edge  of  wings 
is  well  known.  Figures  34  and  35  show  a  comparison  between  the  nonlinear  com¬ 
puted  results  and  the  results  of  the  linearized  panel  method  of  Reference  16. 
Overall  good  correlation  is  achieved  at  these  low  angles  of  attack  except  in  the 
vicinity  of  the  leading  edge.  In  Figure  34  the  linearized  results  depart  sharply 
from  the  nonlinear  at  approximately  95%  span.  The  symmetry  plane  pressures 
correlate  reasonably  well  and  the  loading  or  Cp  is  in  excellent  agreement  over 
most  of  the  wing.  In  Figure  35  the  agreement  in  pressures  on  the  upper  surface 
is  not  as  good  as  for  the  thinner  wing. 

Figures  36  and  37  show  the  same  conditions  as  Figures  34  and  35  but  with 
spanwise  camber.  Linearized  theory  predicts  the  reduction  of  the  upper  surface 
pressures  and  is  in  excellent  agreement  with  the  nonlinear  results.  Some 
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Figure  32  Comparison  of  Lifting  Pressure  Distribution  with  and  without  Camber  for 
14:1  Elliptical  Wing,  2.0,  a  =  3° 


discrepancy  in  the  lower  surface  pressure  exists.  Figure  38  shows  a  compari¬ 
son  of  the  two  methods  at  10"  angle  of  attack  for  a  symmetrical  elliptic  wing. 
Linearized  theory  breaks  down  completely  in  the  region  of  supercritical  flow. 

The  upper  surface  or  leeward  pressures  are  in  fair  agreement  between  the 
shock  and  symmetry  plane.  The  pressures  on  the  windward  side  of  the  wing 
are  underpredictcd  by  the  linearized  method.  The  loading  between  the  shock  and 
the  symmetry  plane  is  also  somewhat  undcrpredicted. 

From  these  comparisons  it  can  be  seen  that  the  breakdown  of  linearized 
methods  is  primarily  localized  to  the  supercritical  regions  of  these  flows.  Lin¬ 
earized  methods  cannot  predict  the  existence  of  a  crossflow  shock  because  this 
effect  is  inherently  nonlinear.  Nonlinear  prediction  techniques  can  be  used  to 
evaluate  the  beneficial  effects  of  aerodynamic  devices  like  spanwise  camber. 

The  lengthy  process  employed  with  the  cone  as  a  starting  solution  had  the 
geometry  first  deform  to  the  desired  conical  cross  section.  The  flowficld  then 
first  converges  on  the  body  solution  and  eventually  on  the  entire  flow  field.  The 
linearized  solution  allows  starting  with  the  desired  conical  cross  section  with 
no  geometrical  deformation.  The  bow  shock  converges  first  and  then  the  entire 
flowfield. 

The  useful  aspect  of  employing  a  linearized  starting  solution  is  that  lower 
free  stream  Mach  numbers  and  larger  delta  angles  can  be  achieved  than  were 
possible  using  the  cone  as  a  starting  solution.  Most  of  the  following  results 
were  obtained  using  a  modified  linearized  potential  flow  solution. 

Figure  39  shows  the  converged  shock  shape,  in  comparison  with  the  Mach 
cone,  for  a  30°  delta  wing  at  MK  =  1 . 8,  a/b  =  22.05.  The  computed  bow  shock 
shows  the  largest  difference  at  the  leading  edge  of  the  wing  where  the  Mach  cone 
is  close  to  the  tip.  Figure  40  shows  the  convergence  history  of  the  bow  shock  and 
body  pressures  for  a  20°  delta  at  M„  =  2.  0.  The  wing  bow  shock  converges  quite 
rapidly  and  the  body  pressure  follows. 
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Using  a  linearized  starting  solution  therefore  seems  plausible,  and  shows 
promise  for  the  computation  of  flow  about  wings.  The  next  step  is  to  generate 
solutions  at  angle  of  attack  using  a  similar  linearized  starting  solution. 

Figure  41a  shows  a  previous  case  recomputed  using  a  modified  linearized 
potential  flow  starting  solution  and  previously  outlined  techniques  to  resolve 
the  vortical  layer  for  the  6:1  elliptic  cone  at  =  1. 97,  a  =  10°.  The  computed 
pressure  distribution,  and  wing  cross  section  are  shown.  Figure  41b  shows  the 
corresponding  crossflow  Mach  number  distribution.  Figure  41c  shows  the  com¬ 
puted  crossflow  streamline  pattern,  along  with  the  sonic  line  and  crossflow  and 
bow  shock  locations. 

Thus  far,  only  delta  wings  with  subsonic  leading  edges  have  been  shown. 
Figure  42  shows  the  pressure  distributions  for  a  20°  delta  wing  with  elliptical 
cross  section  (a/b=  14)  at  a  =  10°  and  Ma0  =  2.  ,  3. ,  4.  and  6.  The  four  Mach 
numbers  correspond  to  a  subsonic  leading  edge,  near-sonic  or  slightly  super¬ 
sonic,  and  supersonic  leading  edges,  respectively.  With  increasing  Mach  num¬ 
ber,  the  crossflow  shock  moves  inboard  toward  the  symmetry  plane.  As  the 
leading  edge  becomes  supersonic,  the  windward  surface  compression  gradient 
becomes  increasingly  steep. . 

Figure  43  shows  the  crossflow  Mach  number  distribution  for  the  same 
cases.  Surprisingly,  the  lowest  Mach  number  shows  the  largest  crossflow  Mach 
number  in  the  vicinity  of  the  leading  edge.  The  crossflow  stagnation  point  moves 
closer  to  the  leading  edge  with  increasing  Mach  number,  causing  less  of  an  ex¬ 
pansion  to  occur  around  the  leading  edge.  Hence,  the  lower  crossflow  Mach 
number  accompanies  increasing  free  stream  Mach  number.  Figure  44  shows 
the  shock  patterns  for  these  Mach  numbers.  At  Mach  2,  a  portion  of  the  leeward 
bow  shock  does  not  show  a  discrete  jump  in  pressure.  With  increasing  Mach 
number,  the  bow  shock  differs  significantly  from  the  Mach  wave,  especially  in 
the  vicinity  of  the  leading  edge. 

The  Mach  2  solution  was  generated  using  a  linearized  potential  flow  start¬ 
ing  solution  and  the  rotated  Mach  wave  as  an  initial  guess  for  the  bow  shock 
shape.  The  Mach  3  solution  was  generated  using  the  Mach  2  Euler  solution  as 
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(b)  CROSSFLOW  MACH  NUMBER  DISTRIBUTION 
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an  initial  guess  for  the  flow  directions  and  the  bow  and  crossflow  shock  shapes. 
The  entire  flowfield  was  scaled  to  the  Mach  3  enthalpy.  The  other  solutions 
were  started  in  a  similar  fashion,  using  the  Mach  3  flowfield.  Figure  45  shows 
the  crossflow  streamline  patterns  for  these  Mach  numbers.  Clearly,  the  flow 
angularity  around  the  leading  edge  in  the  crossflow  plane  is  greatest  for  the 
Mach  2  solution. 

Figure  46a  shows  the  computed  results  at  Mach  2  and  a  =5°  for  three 
planform  angles,  15°,  20°  and,  30°.  Increasing  the  planform  angle  reduces  the 
shock  in  this  case  to  a  compression.  In  general,  for  identical  flow  conditions 
and  similar  geometry,  increasing  the  planform  angle  will  reduce  the  strength 
of  the  shock.  Figure  46b  shows  the  computed  crossflow  Mach  number  distribu¬ 
tions.  The  30°  planform  angle  indicates  a  supercritical  region  without  a  distinct 
crossflow  shock. 

4.  COMPARISON  WITH  NONLINEAR  POTENTIAL  FLOW  SOLUTIONS 

A  new  procedure  was  recently  developed  at  Grumman  to  compute  conical 
supersonic  flows  using  the  complete  nonlinear  potential  flow  equation.  The 
method  was  developed  by  Grossman  (Reference  17).  This  procedure  uses  a 
conical  relaxation  technique  to  obtain  a  solution  of  the  irrotational  Euler  equa¬ 
tions.  With  the  assumptions  of  irrotationality ,  vortical  singularities  and  cross- 
flow  layers  do  not  occur  or  present  any  problems  in  obtaining  converged  solu¬ 
tions.  Hence,  some  of  the  results  of  the  conical  relaxation  method  were  com¬ 
pared  to  the  solutions  obtained  by  integrating  the  rotational  Euler  equations. 

Figures  47  through  51  show  some  of  the  surface  pressure  distribution 
comparisons  for  several  delta  wings  with  elliptical  cross  sections  at  various 
angles  of  attack  and  Mach  numbers.  Figure  47  shows  the  results  for  a  Mach 
1.97,  6:1  ellipse  at  a  =  10°.  The  conical  relaxation  method  captures  the  shocks 
and  the  crossflow  shock  is  usually  smeared  over  two  or  three  grid  points.  The 
comparison  is  quite  good  for  both  the  location  and  pressure  jump  across  the 
shock.  The  conical  relaxation  method  predicts  somewhat  lower  pressures  on 
the  expansion  side  of  the  leading  edge  as  compared  to  the  Euler  solution.  These 
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Figure  48  Comparison  of  Nonlinear  Theories  for  14:1  Ellipse,  M( 
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Figure  51  Comparison  of  Nonlinear  Theories  for  14:1  Ellipse,  Mm 


lower  pressures  result  in  somewhat  higher  crossflow  Mach  numbers  in  the 
vicinity  of  the  leading  edge.  A  slight  discrepancy  in  the  leeward  symmetry 
plane  pressure  is  also  apparent,  while  the  windward  pressure  distribution  is  in 
excellent  agreement. 

Figure  48  shows  a  comparison  for  a  much  thinner  ellipse  (a/b~  14)  at 
Mach  2,  a  =5°.  The  entire  windward  and  the  leeward  portion  of  the  pressure 
distribution  between  the  crossflow  shock  and  symmetry  plane  are  in  excellent 
agreement.  The  conical  relaxation  predicts  a  higher  crossflow  Mach  number 
and  an  expansion  pressure  “spike"  around  the  leading  edge.  The  locations  of 
the  crossflow  shock  exhibit  only  fair  agreement.  Figure  49  shows  the  ,ame 
wing  at  a  =  10°.  Figure  50  shows  an  ellipse  at  a  planform  angle  of  30°  and 
M<o  =  20,  o  =5°.  Excellent  agreement  is  achieved,  with  the  potential  flow  solu¬ 
tion  indicating  a  weak  crossflow  shock. 

Since  the  conical  relaxation  method  is  irrotational,  a  higher  Mach  number 
comparison  was  made  to  determine  if  any  entropy-related  discrepancies  would 
occur.  Figure  51  shows  the  comparison  for  the  delta  wing  of  Figure  49  butat  Mach 
3,  a  =10°.  The  figure  exhibits  the  same  trends  as  the  previous  comparison  with 
no  large  discrepancies  occurring  at  the  higher  Mach  number.  Figures  47  to 
50  correspond  to  subsonic  leading  edges,  Figure  51  corresponds  to  a  slightly 
supersonic  leading  edge. 

The  relaxation  method  predicts  somewhat  higher  Mach  numbers  and  lower 
expansion  pressures  in  the  supersonic  crossflow  region,  but  the  overall  agree¬ 
ment  in  pressure  is  excellent.  The  conical  relaxation  method  adequately  cap¬ 
tures  the  shocks  and  exhibits  all  the  major  nonlinear  characteristics  of  the 
flow.  The  relaxation  method  requires  a  finer  mesh  (60x60)  in  comparison  to 
the  Euler  solutions  (50 x  35),  but  does  not  require  shock  fitting. 

5.  OTHER  CONICAL  SHAPES 

a.  Cones  at  High  Incidence 

Some  examples  of  the  flow  generated  by  circular  cones  at  high  incidence 
were  computed  to  determine  if  significant  lift-off  of  the  vortical  singularity 


could  be  achieved  at  low  supersonic  Mach  numbers.  Figure  52  shows  the  cross- 
flow  streamline  pattern  computed  for  a  5°  cone  semi-apex  angle  at  M„  =  2.0, 
a  =10°.  The  streamlines  clearly  depict  the  lift-off  of  the  vortical  singularity. 
The  circular  cone  streamline  pattern  is  distinct  from  the  elliptic  cone  in  that 
only  one  nodal  singularity  occurs  in  the  leeward  symmetry  plane.  The  stagna¬ 
tion  streamline  wetting  the  body  lies  in  the  windward  symmetry  plane.  Hence, 
the  body  windward  symmetry  plane  point  is  the  saddle  point.  Numerical  diffi¬ 
culties  were  not  encountered  in  computing  this  case  because  of  the  small  cone 
angle.  A  supersonic  crossflow  region  does  not  occur  and  thus  the  absence  of  a 
crossflow  shock.  The  bow  shock  is  weak  and  thus  the  flow  is  not  far  from  ir- 
rotational.  Figure  53  shows  the  result  for  a  12.5°  cone  at  M„=1.8,  a  =25°. 

This  cone  was  computed  in  an  attempt  to  achieve  vortical  lift-off  with  a  cross- 
flow  shock  at  a  lower  Mach  number.  The  streamlines  exhibit  the  normal  coales- 
cent  behavior  at  the  leeward  vortical  singularity  that  is  typical  of  incipient 
lift-off.  Lift-off  was  not  achieved  at  the  angle  of  attack  corresponding  to  twice 
the  cone  angle.  The  crossflow  shock  exhibits  large  curvature  in  its  radial 
shape.  This  may  be  an  extreme  case  of  the  vortical  and  crossflow  layer  prob¬ 
lems  mentioned  earlier.  For  circular  cones,  the  radial  density  of  streamlines 
is  greater  in  the  vicinity  of  the  surface,  thus  indicating  larger  gradients  in 
entropy  and  other  flow  variables.  It  was  necessary  to  use  40  radial  grid  points 
to  have  enough  resolution  to  fit  the  crossflow  shock. 

b.  Wing-Body  Effect 

To  determine  the  interference  effect  of  a  conical  body  on  the  wing  flowfield, 
two  conical  wing-body  geometries  were  computed.  Both  bodies  were  cubically 
faired  to  the  wing  at  a  specified  conical  location  (see  Figure  54).  Both  wing-body 
geometries  were  computed  using  a  converged  Euler  solution  for  the  wing  as  the 
initial  starting  solution.  The  geometry  was  then  deformed  continuously  until  the 
specified  conical  wing-body  geometry  was  achieved,  at  which  point  in  the  march¬ 
ing  technique,  the  wing-body  geometry  was  held  invariant  and  conical  until  con¬ 
vergence  was  achieved.  Figures  55  through  57  show  the  results  for  a  very 


-0.40 


Figure  55  Body  Interference  Effect  on  Surface  Prefture  and  Crossflow 
Shock  for  Well-Blended  Conical  Geometry,  Mx  -  2.0,  a  -  10° 
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Figure  57  Body  Interference  Effect  on  Wing-Body  Crossflow 
Streamline  Pattern  for  Well-Blended  Conical 
Geometry,  =  2.0,  a  =  10° 


blended  wing-body  cross  section  where  the  body  was  faired  to  the  wing  at  a 
position  (denoted  by  SE)  midway  between  the  wing  leading  edge  and  the  symmetry 
plane.  The  body  symmetry  plane  angle  was  4.  5°  (6B)  and  the  extended  center- 
line  angle  of  the  wing  was  1.5°  (6^.  Figure  55  shows  the  resulting  pressure 
distribution  and  crossflow  shock  shape  in  comparison  with  the  converged  wing 
solution  that  was  used  to  start  the  calculation.  The  presence  of  the  faired  body 
pushes  the  location  of  the  crossflow  shock  outboard  toward  the  leading  edge. 

The  crossflow  shock  thus  reacts  earlier  to  the  presence  of  the  low  values  of 
crossflow  that  occur  in  the  leeward  shoulder  region  or  wing-body  juncture. 

The  presence  of  the  body  increases  the  windward  pressures  but  also  somewhat 
increases  the  leeward  pressures  between  the  shock  and  symmetry  plane.  The 
outboard  movement  of  the  shock  reduces  the  lift  in  the  leeward  expansion  re¬ 
gion.  The  crossflow  streamline  pattern  changes  somewhat  in  that  the  leeward 
streamlines  appear  to  coalesce  tangentially  to  the  body  at  the  vortical  singular¬ 
ity,  whereas  the  wing  alone  appeared  to  coalesce  normal  to  the  body. 

Figures  58  through  60  show  the  computed  results  for  a  more  abrupt  (not 
so  blended)  wing-body  combination.  The  body  was  faired  to  the  wing  at  a  5° 
angle  (6E),  with  a  body  centerline  angle  of  5°  (6B)  and  the  wing  thickness  angle 
remaining  at  1.5°.  The  same  general  trends  apply  to  these  results  except  for 
one  interesting  aspect.  The  crossflow  streamline  pattern  (Figure  60)  shows 
that  the  windward  nodal  or  vortical  singularity  has  moved  off  the  symmetry 
plane  and  occurs  at  the  wing  body  juncture.  The  windward  symmetry  plane 
must  remain  a  crossflow  stagnation  point,  but  becomes  a  saddle  point  instead  of 
a  nodal  singularity.  The  crossflow  Mach  number  distribution  confirms  the 
streamline  pattern  in  that  four  instead  of  three  crossflow  stagnation  points 
occur  for  this  geometry.  The  leeward  vortical  singularity  remains  intact,  and 
once  again  the  streamlines  appear  to  converge  tangentially.  The  leeward 
streamlines  thus  locally  behave  in  a  similar  fashion  to  the  leeward  streamlines 
on  a  circular  cone  at  moderate  incidence.  Similar  behavior  of  the  leeward  nodal 
singularity  for  conical  wing -body  configurations  has  been  observed  in  the  non¬ 
linear  potential  flow  solutions,  using  the  methods  developed  in  Reference  17. 
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Figure  58  Body  Interference  Effect  on  Surface  Pressure  and 

Crossflow  Shock  for  Less-Blended  Conical  Geometry, 
-2.0,  a  -  10° 
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Figure  59  Body  Interference  Effect  on  Crossflow  Mach  Number 
Distribution  for  Less-Blended  Conical  Geometry, 

=  2.0,  a  -  10° 
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It  must  be  mentioned  that  the  conformal  mapping  clusters  points  around 
the  leading  edge  of  the  wing.  Hence,  minimal  resolution  was  obtained  in  the 
wing-body  juncture  area  with  50  circumferential  points. 

6.  THREE-DIMENSIONAL  WINGS 

Figure  61  shows  an  isometric  view  of  the  pressure  distributions  com¬ 
puted  for  a  fully  three-dimensional  wing,  along  with  some  representative  span- 
wise  sections  using  a  modified  NASA/Craidon  geometry  program  (Reference 
10).  The  wing  begins  with  a  conical  section  to  Z  =30.  The  starting  solution  used 
was  the  converged  Euler  solution  at  Z  =30  and  was  shown  in  Figure  41.  The 
thickness  was  tapered  parabolically  and  the  camber  was  varied  from  Z  =30  to 
Z  =  150,  resulting  in  a  very  thin  cambered  cross  section  at  the  trailing  edge. 

Figure  62  shows  an  overlay  of  the  same  representative  nondimensional 
pressure  distributions  and  cross  sections.  The  longitudinal  thickness  variation 
expands  the  upper  and  lower  surface  pressures.  The  pressures  in  the  vicinity 
of  the  leading  edge  also  expand.  At  the  thinnest  section,  the  pressure  begins 
to  approach  limiting  pressure  (Cp  =  -  .367  at  M„  =  1.97).  Figure  63  demonstrates 
the  development  of  the  crossflow  Mach  number  in  the  vicinity  of  leading  edge 
for  two  of  the  thinner  cross  sections.  At  Z  =  144.  5,  which  is  near  the  trailing 
edge  of  the  wing,  the  crossflow  approaches  a  Mach  number  of  7.  Certainly, 
this  extreme  expansion  around  the  leading  edge  would  lead  to  boundary  layer 
separation. 
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SECTION  XI 
CONCLUSIONS 

It  is  obvious  from  the  results  presented  here  that  there  are  important 
nonlinear  aspects  of  supersonic  flow  over  wings.  The  results  generated 
during  the  course  of  this  work  prove  that  the  numerical  procedures  used 
here  can  be  successful  in  predicting  these  nonlinear  flowfields.  There¬ 
fore  the  code  which  has  been  developed  can  be  a  valuable  aid  in  the  design 
of  supersonic  wings  and  wing/body  combinations. 

The  computational  time  (code  running  time)  required  to  generate 
these  flowfields  seems  to  make  the  code  a  difficult  design  tool  to  use. 

From  the  results  generated  thus  far,  it  seems  that  a  minimum  of  about  one 
hour  C.P.U.  time  is  required  to  predict  the  flowfield  about  a  typical 
supersonic  wing. 

This  number  is  unacceptable  if  the  design  process  is  to  cycle  through 
many  wings  to  achieve  certain  design  goals.  A  significant  reduction  of 
this  running  time  will  require  significant  changes  in  the  computational 
technique  used  in  this  work. 

During  the  course  of  this  work  an  alternative  approach  has  been  con¬ 
sidered.  Specifically,  the  wing  designer  would  specify  the  design  goals 
and  a  computational  procedure  would  develop  a  wing  geometry  to  meet 
these  goals.  Computational  procedures  which  solve  this  inverse  problem 
have  been  reported  quite  extensively  over  the  past  few  years.  Most  of 
these  schemes  involve  iterating  the  wing  geometry,  for  example,  to  achieve 
certain  design  goals.  This  procedure  would  necessarily  involve  a  solution 
to  the  direct  problem  which,  in  this  case,  would  be  unacceptable  because 
of  the  lengthy  computational  times  involved  in  the  direct  problem. 

For  the  simplest  of  the  inverse  problems  (i.e.  specify  the  surface 
pressure  distribution  and  generate  a  wing  geometry),  a  scheme  has  been 
studied  during  the  course  of  this  work  which  could  have  the  same  running 
time  as  the  direct  problem  (1  hr.  C.P.U. ).  It  would  involve  inverting 
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the  surface  pressure  equation  (equation  42)  to  compute  the  surface  geom¬ 
etry  which  corresponds  to  a  given  pressure  distribution.  This  is  a  direct 
scheme  for  computing  the  inverse  (i.e.  design)  problem  and  therefore  the 
computational  time  would  be  reduced  significantly  over  an  iterative  tech- 


APPENDIX  -  COEFFICIENTS  OF  TRANSFORMED  EULER  EQUATIONS 


A,  =w2  -yT  =w2  -a2 


A2=ez  +  (exu  +  0yv)  w/A, 

A  s  =  dz+Ox\i/w  +  dyv/w 

A4  =  (0xu  +  0yv)/Aj 

A5=02  +  (rxu+ryv)w/A, 

A6  »  r2  +  rx  u/w  +  ryv/w 

A7  =  (rxu+ryv)/At 

all  =Xj  +XrA5  +X^A2 

a12=Xxrw/At 

a13=Xyrw/A1 

ai4  =  —  YCXjA  7  +  X$A  4) 

bii^Yj  +  Y^  +  Y.Aj 

b12  =  ywYx/Ai 

bi3 =ywYy/A1 

l>i4  =  -T(YrA7  +  YeA4) 

a21=XxT/w 

a22  =X.  +  XrAg  +  X*A3 
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Figure  62  Pressure  Distributions  and  Crossflow  Shocks  at  Various  Stations 
of  Fully  3  D  Wing,  -  1 .97,  o  ■=  10" 


